The following is the presentation that I gave at the Computational in Economics and Finance (CEF) conference 2014 in Oslo.

Note that all Julia packages used in this notebook are installable directly with the package manager. The only exception is the `splines`

package which is installed with: `Pkg.clone("https://github.com/EconForge/splines.git")`

.

# Overview of Julia¶

Julia is a new programming language for technical computing, created at the MIT in 2012. Its main features are:

- Syntax close to MATLAB/Octave, with ideas also borrowed from Python, Ruby, Lisp
- Good performance, close to compiled languages, thanks to just-in-time (JIT) compiler
- Free software / open source
- Designed for parallel and distributed computation
- Easily interfaceable with other languages (C, Fortran, Python)
- Multiple dispatch: ability to define function behavior across many combinations of argument types
- Large number of package extensions, lively community

# Example 1: Hodrick-Prescott filter¶

Here we show how to implement a Hodrick-Prescott filter in Julia. Recall that, given a timeseries $y_t$ and a smoothing parameter $\lambda$, the filtered trend $\tau_t$ is defined by: $$\min_{\tau}\left(\sum_{t = 1}^T {(y_t - \tau _t )^2 } + \lambda \sum_{t = 2}^{T - 1} {[(\tau _{t+1} - \tau _t) - (\tau _t - \tau _{t - 1} )]^2 }\right)$$

The first-order conditions of this minimization problem can be expressed by the following linear system: $$\begin{pmatrix} 1+\lambda & -2\lambda & \lambda & & & & \\ -2\lambda & 1+5\lambda & -4\lambda & \ddots & & & \\ \lambda & -4\lambda & 1+6\lambda & \ddots & \ddots & & \\ & \ddots & \ddots & \ddots & \ddots & \ddots & \\ & & \ddots & \ddots & 1+6\lambda & -4\lambda & \lambda \\ & & & \ddots & -4\lambda & 1+5\lambda & -2\lambda \\ & & & & \lambda & -2\lambda & 1+\lambda \end{pmatrix} \begin{pmatrix} \tau_1 \\ \vdots \\ \vdots \\ \tau_t \\ \vdots \\ \vdots \\ \tau_T \end{pmatrix} = \begin{pmatrix} y_1 \\ \vdots \\ \vdots \\ y_t \\ \vdots \\ \vdots \\ y_T \end{pmatrix} $$

The following function implements the Hodrick-Prescott filter by computing the filtered trend, given a timeseries and a smoothing parameter:

```
function hp_filter(y::Vector{Float64}, lambda::Float64)
n = length(y)
@assert n >= 4
diag2 = lambda*ones(n-2)
diag1 = [ -2lambda; -4lambda*ones(n-3); -2lambda ]
diag0 = [ 1+lambda; 1+5lambda; (1+6lambda)*ones(n-4); 1+5lambda; 1+lambda ]
D = spdiagm((diag2, diag1, diag0, diag1, diag2), (-2,-1,0,1,2))
D\y
end
```

The `spdiagm`

function is used to create a matrix by specifying its nonzero diagonals. The right-divide operator (antislash) is used to compute the solution of the linear system.

Features to be noted:

- Matrices can be easily constructed with brackets
- Function arguments are typed: this is optional most of the time, but improves performance
- The type
`Float64`

represents double precision floating point numbers (they occupy 64 bits in memory). It is the default type when manipulating real numbers. - There is no output argument (it is the last computation that is returned as output)
- The implicit multiplication operator in monomials (
`2lambda`

) - One can construct tuples using a comma-separated list within parenthesis (see the 5-tuples as first and second argument of
`spdiagm`

).

Now we construct a timeseries as the simulation over 40 periods of a random walk with gaussian innovation:

```
y = cumsum(randn(40))
```

And we extract the filtered trend, with a smoothing parameter of 1600 (typical for quarterly frequency):

```
tau = hp_filter(y, 1600.)
```

Note that there is a dot at the end of `1600`

, in order to emphasize the fact that it is a floating point. If the dot is omitted, the number is interpreted as an integer, and this leads to an error. This is a consequence of the typing system of Julia, we'll return to that later.

```
hp_filter(y,1600)
```

Now we want to plot the result, using the `PyPlot`

module which is based on Python's *matplotlib*. Note that a module needs to be loaded only once per Julia session.

```
using PyPlot
```

Finally let's do the plot:

```
plot(1:40, y, "r", 1:40, tau, "b")
title("HP-filtered random walk")
```

# Example 2: yield-to-maturity computation¶

Suppose that we want to compute the yield-to-maturity (YTM) of a bond with pays a yearly coupon $C$. The principal is 100, which is due at date $T$ (expressed in years, today is date 0). The bond is priced today at $P$ (excluding accrued interest).

First consider the case when $T$ is an integer. For a given annual discount rate $\delta$, the net present value of the flow of payments is the following: $$NPV(\delta) = \frac{100}{(1+\delta)^T}+ \sum_{t=1}^{T} \frac{C}{(1+\delta)^t}$$

If $T$ is not an integer, and if $\lfloor T \rfloor$ is the biggest integer less than $T$, then the NPV is: $$NPV(\delta) = \frac{1}{(1+\delta)^{T-\lfloor T \rfloor}} \left( \frac{100}{(1+\delta)^{\lfloor T \rfloor}}+ \sum_{t=1}^{\lfloor T \rfloor} \frac{C}{(1+\delta)^t}\right)$$

The YTM is the value $\delta^*$ that equalizes the NPV with the bond price (including accrued interests), *i.e.* such that:
$$P + (1-T+\lfloor T \rfloor)C = NPV(\delta^*)$$

We now implement this computation in Julia. First, we load two packages that we will need for our purpose. The first is `DateTime`

, useful for time computations, and the second is `Roots`

, which provides a nonlinear zero finder for univariate functions via the `fzero`

function.

```
using Datetime
using Roots
```

Now the function with computes the YTM, in percentage points:

```
function ytm(price::Float64, maturity::Date, coupon::Float64)
function npv_minus_price(yield)
# Compute the NPV sum by backward induction
npv = 100+coupon
d = maturity - today()
while d >= year(1)
npv = coupon + npv/(1+yield/100)
d -= year(1)
end
# Adjustment if the maturity horizon is not an integer number of years
npv /= (1+yield/100)^(integer(d)/365)
# Difference between NPV and market price (with accrued interest)
npv - (price + (1-integer(d)/365)*coupon)
end
fzero(npv_minus_price, coupon)
end
```

The nested function `npv_minus_price`

computes the difference between the NPV for a given yield, and the market price including accrued interest. Then the `fzero`

function computes a zero of the `npv_minus_price`

function, which is the YTM. The zero solver needs an initial value for starting the algorithm, and we use the coupon rate for that purpose.

Features to be noted:

- Functions can be nested. Inside functions inherit from variables of the enclosing function
- The
`while`

/`end`

control structure, similar to MATLAB's - Updating arithmetic operators (
`-=`

,`/=`

), which do both a computation and a variable assignment (*e.g.*`x /= 2`

is equivalent to`x = x/2`

) - Type converters:
`year()`

,`integer()`

- Fields of a structure can be accessed with a dot syntax, as in MATLAB (
*e.g.*`r.zero`

)

```
ytm(101.4, date("2015-09-01"), 2.)
```

The YTM of a bond maturing on September 1st 2015, with a 2% coupon rate, and currently priced 101.4, is 0.82%.

# Types¶

## Concrete bits types¶

The default integer type is 32-bit or 64-bit, depending on your machine hardware:

```
typeof(1)
```

`Int`

is an alias for the default integer type:

```
Int
```

The default floating point type is 64-bit (equivalent of Fortran's `double precision`

or C's `double`

):

```
typeof(1.0)
```

## Abstract types¶

Abstract types do not correspond to actual objects in memory. They are generic categories that include (generally several) concrete types. For exemple, the abstract type `Integer`

includes both 32-bit and 64-bit concrete integer types (one says that `Int32`

and `Int64`

are subtypes of `Integer`

):

```
Int64 <: Integer
```

```
Int32 <: Integer
```

But 64-bit floating points are not a subtype of `Integer`

:

```
Float64 <: Integer
```

## Tuples¶

Tuples are made of several values enclosed within parentheses and separated by commas. They are useful when you want a function to return several values (possibly of different types).

```
tpl = (1,true,"foo")
```

```
typeof(tpl)
```

Tuple elements are indexed with brackets, and indices start at 1:

```
tpl[2]
```

Tuples can be used for multiple variable assignments:

```
(a, b, c) = tpl
```

A special case is the tuple containing only one value. Note the syntax with parenthesis and a trailing comma.

```
typeof((1,))
```

## Array types, matrices¶

A vector of floats:

```
[1., 2., 3.]
```

A vector of integers:

```
[1, 2, 3]
```

Note that `Array{T,n}`

is a parametric type: the parameter `T`

determines the type of objects inside the array, and the parameter `n`

determines the dimensionality.

`Vector{T}`

is an alias for `Array{T,1}`

.

A matrix of floats:

```
[1. 2.; 3. 4.]
```

`Matrix{T}`

is an alias for `Array{T,2}`

.

An unitialized array is created in the following way (beware, elements are not guaranteed to be zero, and can be anything):

```
X = Array(Float64, 10, 5)
```

The size of the array is returned as a tuple:

```
size(X)
```

Array elements are indexed with brackets, and indices start at 1:

```
X[10,2]
```

Multiplication and division operators on matrices perform linear algebra operations:

```
A = [ 1. 2.; 3. 4. ]
B = [ 5. 6. ]'
A * B
```

```
A \ B
```

Element-by-element operations are available under `.+`

, `.-`

, `.*`

, `./`

:

```
A .* A
```

It is also possible to create vectors using *comprehension* syntax, which is similar to the mathematical syntax for describing sets:

```
[ 2x^2 + 1 for x = 1:5 ]
```

## Composite types¶

Composite types are records (also called structures or objects in other languages) containing several named fields. For example, one can create a `Timeseries`

type, which is the aggregation of a start date, a frequency, and an array of data:

```
type Timeseries
start::Date
frequency::Integer
data::Vector{Float64}
end
```

```
x = Timeseries(date("2014-06-21"), 12, [ 1., 2., 3. ])
```

It is possible to read and write members of the composite type using the dot operator:

```
x.data
```

```
x.frequency = 4
```

It is also possible to create composite types whose members cannot be modified, by replacing the keyword `type`

with `immutable`

. Such types are generally more efficient, especially if they are packed inside arrays.

## Type annotations¶

Type annotations can be used at different places:

- in the declaration of composite types (see above)
- for qualifying function arguments (see above)
- to assert that something is of a given type
- to indicate the type of a variable (only local variables in functions, not for global variables)

Note that type annotations are always optional (except for function arguments when doing multiple dispatch, see below). Type annotations can increase performance when used judiciously, see below.

Here is an example of type assertions:

```
1::Integer
```

```
1.0::Integer
```

Here is an example of typed local variables:

```
function f(x::Float64)
y::Integer = floor(x)
2y
end
```

```
f(1.)
```

# Functions, methods¶

## Functions¶

Functions can be constructed using three syntaxes:

- with the
`function`

keyword - using a compact “assignement form” when the function body is simple
- as anonymous function (
*i.e.*not associated to a function name), with the arrow syntax (`->`

)

Here is an example of compact function declaration:

```
g(x) = 2x^2-3x+1
```

```
g(5)
```

Here is an example of anonymous polynomial function, for which we search a zero:

```
fzero(x -> 2x^2-3x, 1.)
```

Functions are also objects of the language, and have a type. This can be checked on an anonymous function (*i.e.* a function which is not given a name):

```
typeof(x -> x+1)
```

One interesting feature when declaring functions is the ability to declare optional arguments, that will be referenced by a keyword (like in R). Those optional arguments have to be listed after a semicolon in the function declaration, and should be given a default value:

```
function utility(c; gamma = 2.0)
c^(1-gamma)/(1-gamma)
end
utility(2)
```

```
utility(2, gamma = 3.0)
```

## Methods and multiple dispatch¶

In Julia, a given function can behave differently depending on the types of its arguments. In this case, the different implementations are called *methods*, and the process by which a given method is chosen depending on the argument types is called *multiple dispatch*. For example, the square root function has several methods:

```
methods(sqrt)
```

Now we define addition over two timeseries, using the `Timeseries`

type defined above. For the sake of simplicity, we require that the two timeseries have the same start date, the same frequence and the same length.

```
function +(a::Timeseries, b::Timeseries)
if a.start != b.start
error("Arguments do not have the same start date")
end
if a.frequency != b.frequency
error("Arguments do not have the same frequency")
end
if length(a.data) != length(b.data)
error("Arguments do not have the same data length")
end
Timeseries(a.start, a.frequency, a.data .+ b.data)
end
```

Now we can add two timeseries using the most natural syntax:

```
Timeseries(today(), 12, [1., 2., 3]) + Timeseries(today(), 12, [4., 6., 9.])
```

Note that multiple dispatch works not only on existing operators, but also on user-defined functions.

# Performance¶

Julia has a good performance, close to compiled languages like Fortran or C, and often much better than MATLAB. This is evidenced by the micro benchmarks performed on several simple algorithms accross several languages (including C, Fortran, Python, MATLAB, R).

The good performance of Julia comes from its use of a just-in-time (JIT) compiler: all the code that you write in Julia is invisibly compiled on-the-fly by the interpreter.

To be effective, the JIT compiler needs to determine types of variables and functions arguments. Most of the time, types are automatically inferred by Julia in an efficient way. But sometimes the compiler needs some help, and this is where type annotations become useful when put at strategic place. A typical Julia workflow is therefore to prototype a first version of the code without type annotations; then, during a second pass, the code is optimized by putting type annotations at strategic places.

As a general rule, type annotations are useful within composite type declarations (*i.e.* within `type`

and `immutable`

blocks), in function arguments, and in some function local variables. Global variables (*i.e.* variables declared/assigned outside functions) should generally be avoided, or should be declared as constants (using the `const`

keyword, like in `const a = 1`

), because they currently cannot be typed.

For more details on optimization, the Julia manual has some performance tips.

# Where to get more information¶

Here are some links where you can find documentation or extra information about Julia:

- Project homepage
- Julia documentation: an overview of all features of the language and of all primitive functions in the base library
- External packages: list of all add-on package

# Exercise: Solving a rational-expectations model with value function iteration¶

We are going to solve a model by Aguiar and Gopinath (2006): “Defaultable debt, interest rates and the current account”, *Journal of International Economics*, 69(1). It is model II of the paper, which is a model of strategic sovereign default with a stochastic growth trend. The solution method is value function iteration, with approximation of value function with cubic splines.

## Preliminary: Gauss-Hermite quadrature¶

For solving the model, we will need to integrate expectations with respect to Gaussian shocks. We will therefore use the Gauss-Hermite quadrature, which is well-suited to that case. The quadrature consists in approximating an integral by a finite weighted sum. More precisely, the Gauss-Hermite quadrature is: $$\int_{-\infty}^{+\infty} f(x)e^{-x^2} dx \approx \sum_{i=1}^{n} w_i f(x_i)$$ where $n$ determines the precision of the approximation (the higher $n$, the more precise), the $x_i$ are the roots of the $n$-th order Hermite polynomial $H_n(x)$, and the associated weights are: $$w_i = \frac{2^{n-1} n! \sqrt{\pi}}{n^2 [H_{n-1}(x_i)]^2}$$

In turn, the Hermite polynomials are defined by the following recursive equations (where $H'$ denotes the derivative of $H$): $$H_0(x) = 1$$ $$H_n(x) = 2xH_{n-1}(x) - H'_{n-1}(x)$$

We are now going to implement the computation of the points and weights as a function of $n$. For that purpose, we will use the `Polynomial`

package, whose documentation can be seen there.

First, load the package.

```
using Polynomial
```

Now, write a function that, given $n$, creates $H_n$. You will need the `Poly`

function to create basic polynomials, the `polyder`

function to compute the derivative of a polynomial. Standard addition and multiplication work on polynomials as expected. Use a `for`

loop to do the recursion.

```
function compute_hermite_polynomial(n)
P = Poly([1])
const x = Poly([1; 0])
for i = 1:n
P = 2x*P - polyder(P)
end
P
end
```

Verify that $H_5(x) = 32x^5-160x^3+120x$:

```
compute_hermite_polynomial(5)
```

Now, write a function that, given $n$, returns a tuple consisting of the two vectors $x_i$ and $w_i$. Use the `roots`

function to compute the roots of a polynomial, the `polyval`

function to evaluate a polynomial, and the `factorial`

function. You can use a comprehension to create the vector of weights. Also, you can use a type annotation to indicate to Julia that the roots of the polynomial are real, and not complex.

```
function compute_hermite_points(n)
P = compute_hermite_polynomial(n)
P2 = compute_hermite_polynomial(n-1)
points::Vector{Float64} = roots(P) # We know that the roots are real
weights = [ 2^(n-1)*factorial(n)*sqrt(pi)/n^2/polyval(P2, points[i])^2 for i = 1:n ]
points, weights
end
```

Now, verify that, for $n=3$, one has: $(x_1,w_1) \approx (-1.22,0.29)$, $(x_2,w_2) \approx (0,1.18)$, $(x_3,w_3) \approx (1.22,0.29)$:

```
compute_hermite_points(3)
```

## The model¶

We consider a country with a representative agent, producing and consuming a single homogeneous good. The production stream is exogenous and has a stochastic growth trend.

At every period, the difference between production and consumption is financed on international markets. The country therefore accumulates debt (or assets) which is short-term (one period), and can decide to default or repay at every period. If it repays, access to financial markets continues. If it defaults, it is barred from financial markets (and therefore, being in autarky, has to consume what it produces), and moreover suffers from a penalty (proportional to GDP); exclusion from financial markets can last several periods, and terminates with a given probability at each period.

International investors are risk-neutral and in perfect competition. The interest rate that the country faces is therefore equal to the international riskless rate (assumed to be constant), plus a risk premium equal to the model-consistent risk of default in next period.

Let's now describe the model equations. Since the model has a stochastic trend, equations and variables are in detrended form. The two state variables are $y_t$ the (log of the gross) growth rate and $a_t$ the detrended assets due at date $t$.

The law of motion of $y_t$ is: $$y_t = (1-\rho_g)\mu_y +\rho_g y_{t-1} + \varepsilon_t$$ where $\varepsilon_t \leadsto \mathcal{N}(0, \sigma_g)$ and $\mu_y = ln(\mu_g) - \frac{\sigma_g^2}{2(1-\rho_g^2)}$.

Detrended GDP as a function of $y_t$ is: $$Q_t = \frac{e^{y_t}}{\mu_g}$$

We now describe the value function of the country. The instantaneous utility function is: $$u(c) = \frac{c^{1-\gamma}}{1-\gamma}$$ and the discount factor is $\beta$.

$V^G(a_t, y_t)$ is the value function if the country repays, $V^B(y_t)$ is the value function under autarky (after a default), and $V(a_t,y_t)$ is the optimal choice between repayment and default, which is given by: $$V(a_t,y_t) = \max\{ V^G(a_t,y_t), V^B(y_t) \}$$

The value function under autarky depends on $\delta$ (the penalty) and $\lambda$ (the probability of redemption): $$V^B(y_t) = u\left[(1-\delta)Q_t\right] + \beta e^{y_t(1-\gamma)}\left[\lambda \mathbb{E}_t V^G(0, y_{t+1}) + (1-\lambda) \mathbb{E}_t V^B(y_{t+1})\right]$$

The value function under repayment depends of $q(y_t, a_{t+1})$ which is the price at which investors accept to buy bonds due tomorrow: $$V^G(a_t, y_t) = \max_{a_{t+1}} \left\{u\left[Q_t + a_t - a_{t+1} q(y_t, a_{t+1}) e^{y_t}\right] + \beta e^{y_t(1-\gamma)} \mathbb{E}_t V(a_{t+1}, y_{t+1})\right\}$$

Finally, the price function of bonds is: $$q(y_t, a_{t+1}) = \frac{\mathbb{P}(y_{t+1} \geq y^*(a_{t+1}) | y_t)}{1+r}$$ where $r$ is the riskless interest rate, and $y^*(a_{t+1})$ is the default threshold defined by $V^G(a_{t+1}, y^*(a_{t+1})) = V^B(y^*(a_{t+1}))$.

## Calibration¶

We are going to use the following (quarterly) calibration: $\gamma = 2$, $r = 0.01$, $\beta=0.8$, $\lambda=0.1$, $\delta=0.02$, $\mu_g = 1.006$, $\sigma_g = 0.03$, $\rho_g = 0.17$.

Assign each parameter to a constant variable of the same. Also create a variable for $\mu_y$, using the relationship shown above.

```
const gam = 2.0
const r = 0.01
const beta = 0.8
const lambda = 0.1
const delta = 0.02
const mu_g = 1.006
const sigma_g = 0.03
const rho_g = 0.17
const mu_y = log(mu_g)-0.5*sigma_g^2/(1-rho_g^2)
```

## Utility and GDP¶

Using the short form syntax, create functions for computing the instantaneous utility (as a function of consumption) and GDP (as a function of $y_t$):

```
util(c) = c^(1-gam)/(1-gam)
```

```
gdp(y) = exp(y)/mu_g
```

## Grid¶

As explained above, the two state variables are $y_t$ and $a_t$. Since we are going to interpolate with cubic splines, we have to decide of a grid on which to make the approximation. We will use the interval $[-0.3,0]$ for $a_t$, with 15 points. And the interval $\left[\mu_y-\frac{6\sigma_g}{\sqrt{1-\rho_g^2}}, \mu_y+\frac{6\sigma_g}{\sqrt{1-\rho_g^2}}\right]$ for $y_t$, with 10 points.

Create variables `min_a`

, `max_a`

, `ngrid_a`

and `min_y`

, `max_y`

, `ngrid_y`

. Then create vectors `grid_a`

and `grid_y`

, using the `linspace`

function.

```
const min_a = -0.3
const max_a = 0.0
const ngrid_a = 15
const min_y = mu_y - 6*sigma_g/sqrt(1-rho_g^2)
const max_y = mu_y + 6*sigma_g/sqrt(1-rho_g^2)
const ngrid_y = 10
const grid_a = linspace(min_a, max_a, ngrid_a)
const grid_y = linspace(min_y, max_y, ngrid_y)
```

Now, let's give initial values for $V^G$ and $V^B$ on the grid points. We will use $V^B(y_t) = u((1-\delta)Q_t)$ and $V^G(a_t, y_t) = u(Q_t + a_t)$.

Create a vector `VB_points`

and a matrix `VG_points`

initialized with these values. For `VB_points`

, you can use a comprehension. Declare both variables as constant; note that this means that the variable will always point to the same array, but this does not prevent modification of array elements.

```
const VB_points = [ util((1-delta)*gdp(y)) for y in grid_y ]
const VG_points = Array(Float64, ngrid_a, ngrid_y)
for iy = 1:ngrid_y, ia = 1:ngrid_a
VG_points[ia, iy] = util(gdp(grid_y[iy])+grid_a[ia])
end
```

## Splines¶

To interpolate the value function between grid points, we will use cubic splines, available in the `splines`

package.

```
using splines
```

The following example shows how to approximate the absolute function by a 5-points spline over $[-1, 1]$. For 2-dimensional approximation, the principle is the same: arguments to `filter_coeffs`

and `eval_UC_spline`

should be vectors of 2 elements (or matrices for the point values or the coefficients).

```
# Create a small grid
mygrid = linspace(-1., 1., 5)
# Value of the function at grid points
mypoints = Float64[ abs(x) for x in mygrid ]
# Compute coefficients of the spline interpolation
mycoeffs = filter_coeffs([-1.], [1.], [5], mypoints)
# Create a function that computes the interpolation at any point, using the coefficients
interp(x) = getindex(eval_UC_spline([-1.], [1.], [5], mycoeffs, [x]), 1)
# Plot the absolute function and its spline approximation on a finer grid
myfinergrid = linspace(-1., 1., 100)
plot(myfinergrid, [ interp(x) for x in myfinergrid ], "r", myfinergrid, abs(myfinergrid), "b")
```

Now, create constant arrays `VG_coeffs`

and `VB_coeffs`

for storing the splines approximations, using the `filter_coeffs`

function and the arrays `VG_points`

and `VB_points`

.

```
const VG_coeffs = filter_coeffs([min_a, min_y], [max_a, max_y], [ngrid_a, ngrid_y], VG_points)
const VB_coeffs = filter_coeffs([min_y], [max_y], [ngrid_y], VB_points)
```

Write functions `VG(a, y)`

and `VB(y)`

that evaluate the spline at any point in the state space. *Warning*: for `VG`

, the last argument in the call to `eval_UC_spline`

must be `[a, y]'`

(note the transpose, needed to create a $1\times 2$ matrix).

```
VG(a, y) = getindex(eval_UC_spline([min_a, min_y], [max_a, max_y], [ngrid_a, ngrid_y], VG_coeffs, [a, y]'), 1)
```

```
VB(y) = getindex(eval_UC_spline([min_y], [max_y], [ngrid_y], VB_coeffs, [ y ]), 1)
```

Verify that `VB(0.01)`

is approximatively `-1.01631`

and that `VG(-0.01, 0.02)`

is approximatively `-0.9959`

.

```
VB(0.01)
```

```
VG(-0.01, 0.02)
```

Finally, write the function `V`

that is the maximum of `VB`

and `VG`

:

```
V(a, y) = max(VG(a, y), VB(y))
```

## Bond price function¶

First, here is an example of the probability that a variable following $\mathcal{N}(\mu_y, \sigma_g)$ is less than $2\mu_y$, using the `Distributions`

package:

```
using Distributions
cdf(Normal(mu_y, sigma_g), 2mu_y)
```

Now write the bond price function `q(y, a1)`

(where `a1`

denotes $a_{t+1}$). The idea is to compute tomorrow's default threshold $y^*$ (using `fzero`

on the difference between `VG`

and `VB`

), and then to compute the corresponding default probability with the `cdf`

and `Normal`

functions.

In order to speed up the computation, you can treat special cases separately:

- if the assets are positive, the risk premium is zero
- if for a large value of $y_{t+1}$ (like 6 standard deviations above the mean), default is still preferable to repayment, then the bond price is approximatively zero
- if for a small value of $y_{t+1}$ (like 6 standard deviations below the mean), repayment is still preferable to default, then the risk premium is approximatively zero

```
function q(y, a1)
if a1 >= 0
return 1/(1+r)
end
const cond_mean = (1-rho_g)*mu_y + rho_g*y
const width = 6.0
const upper_y1 = cond_mean + width*sigma_g
const lower_y1 = cond_mean - width*sigma_g
if VG(a1, lower_y1) >= VB(lower_y1)
return 1/(1+r)
end
if VG(a1, upper_y1) <= VB(upper_y1)
return 0.0
end
y1_threshold = fzero(y1 -> VG(a1, y1) - VB(y1), lower_y1, upper_y1)
Erepay1 = 1-cdf(Normal(cond_mean, sigma_g), y1_threshold)
Erepay1/(1+r)
end
```

## Bellman equation¶

We are now going to precompute points and weights for quadrature. Define a constant containing the number of quadrature points, which we will set to 16. Then, compute the hermite points and weights using the function created previously.

```
const gauss_npoints = 16
gauss_points, gauss_weights = compute_hermite_points(gauss_npoints)
```

Now, if $x_i$ and $w_i$ are the hermite points, and given that $\varepsilon_t \leadsto \mathcal{N}(0, \sigma_g)$, the following approximation holds (using the change of variable $y = \sigma_g\sqrt{2}x$): $$\mathbb{E}_t f(\varepsilon_{t+1}) = \int_{-\infty}^{+\infty} \frac{1}{\sigma_g\sqrt{2\pi}}f(y)e^{-\frac{y^2}{2\sigma_g^2}}dy = \int_{-\infty}^{+\infty} \frac{1}{\sqrt{\pi}} f(\sigma_g\sqrt{2}x)e^{-x^2}dx \approx \sum_{i=1}^n \frac{w_i}{\sqrt{\pi}} f(\sigma_g\sqrt{2} x_i) = \sum_{i=1}^n w'_i f(x'_i)$$ where $w'_i=\frac{w_i}{\sqrt{\pi}}$ and $x'_i=\sigma_g\sqrt{2} x_i$.

Compute the vector of $x'_i$ and $w'_i$.

```
gauss_weights /= sqrt(pi)
gauss_points *= sqrt(2) * sigma_g
```

Write the function `VB_next(y)`

that, given $y_t$, computes $u\left[(1-\delta)Q_t\right] + \beta e^{y_t(1-\gamma)}\left[\lambda \mathbb{E}_t V^G(0, y_{t+1}) + (1-\lambda) \mathbb{E}_t V^B(y_{t+1})\right]$ (the right hand side in the Bellman equation defining $V^B$).

```
function VB_next(y)
const cond_mean = (1-rho_g)*mu_y + rho_g*y
EVB1 = 0.0
EVG1 = 0.0
for i = 1:gauss_npoints
EVB1 += gauss_weights[i]*VB(cond_mean + gauss_points[i])
EVG1 += gauss_weights[i]*VG(0.0, cond_mean + gauss_points[i])
end
util((1-delta)*gdp(y))+beta*exp(y*(1-gam))*(lambda*EVG1 + (1-lambda)*EVB1)
end
```

Write `VG_objective(a1, a, y)`

which, given $a_{t+1}$, $a_t$ and $y_t$, computes $u\left[Q_t + a_t - a_{t+1} q(y_t, a_{t+1}) e^{y_t}\right] + \beta e^{y_t(1-\gamma)} \mathbb{E}_t V(a_{t+1}, y_{t+1})$ (the maximization objective in the right hand side of the Bellman equation defining $V^G$).

```
function VG_objective(a1, a, y)
const cond_mean = (1-rho_g)*mu_y + rho_g*y
EV1 = 0
for i = 1:gauss_npoints
EV1 += gauss_weights[i]*V(a1, cond_mean + gauss_points[i])
end
util(gdp(y)+a-exp(y)*a1*q(y,a1)) + beta*exp(y*(1-gam))*EV1
end
```

## Main loop¶

First, load the `Optim`

package that will be used for optimization of the value function objective (see the package documentation).

```
using Optim
```

Here is a simple example of minimization:

```
optres = optimize(x -> 2x^2+1, -2., 2.)
optres.f_minimum
```

Create constant variables to store the maximum number of iterations (set to 1000), and the convergence criterion expressed as a tolerance (set to $10^{-6}$).

Then create an iteration counter (initialized to 1), and a variable to store the current error (initialized to $+\infty$).

```
const max_it = 1000
const convergence_tol = 1e-6
it = 1
err = Inf
```

Create two uninitialized arrays `VG_points_next`

and `VB_points_next`

, respectively of the same size than `VG_points`

and `VB_points`

, that will be used to temporarily store the values of value functions at grid points for the next iteration.

```
const VG_points_next = Array(Float64, ngrid_a, ngrid_y)
const VB_points_next = Array(Float64, ngrid_y)
```

Now create the main computational loop:

- use a
`while`

construct, which will iterate as long as the error is greater than the tolerance, and the iteration counter less than the maximum iterations - in an embedded loop over the $y_t$ grid:
- compute the value of $V^B$ at next iteration using
`VB_next`

, and store it in`VB_points_next`

- in an embedded loop over the $a_t$ grid, optimize
`VG_objective`

over $a_{t+1}$, and store the result in`VG_points_next`

- compute the value of $V^B$ at next iteration using
- update the error, which is the max of the (infinite norm) distance between
`VG_points`

and`VG_points_next`

on one hand, and`VB_points`

and`VB_points_next`

on the other hand - using
`copy!`

, update`VG_points`

and`VB_points`

with the new values - using
`copy!`

and`filter_coeffs`

, update the spline coefficients - using
`print`

, display the iteration counter and the current error - increment the iteration counter

```
while err > convergence_tol && it <= max_it
for iy = 1:ngrid_y
const y = grid_y[iy]
VB_points_next[iy] = VB_next(y)
for ia = 1:ngrid_a
a = grid_a[ia]
optres = optimize(a1 -> -VG_objective(a1, a, y), min_a, max_a)
VG_points_next[ia,iy] = -optres.f_minimum
end
end
# Compute error
err = max(maximum(abs(VB_points_next-VB_points)),
maximum(abs(VG_points_next-VG_points)))
# Copy
copy!(VB_points, VB_points_next)
copy!(VG_points, VG_points_next)
copy!(VG_coeffs, filter_coeffs([min_a, min_y], [max_a, max_y], [ngrid_a, ngrid_y], VG_points))
copy!(VB_coeffs, filter_coeffs([min_y], [max_y], [ngrid_y], VB_points))
# Output
print("Iteration ", it, ": error=", err, "\n")
it += 1
if it == max_it+1
warn("Maximum iterations reached")
end
end
```

Verify that for $y_t=0$, the default threshold is 21% of debt to (quarterly) GDP:

```
fzero(a -> VG(a, 0.) - VB(0.), 0.)/gdp(0.)
```

Copyright © 2014 Sébastien Villemot `<sebastien.villemot@sciencespo.fr>`

License: Creative Commons Attribution-Share Alike 4.0 or GNU General Public License version 3 or later

Go Top