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:

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:

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

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.

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.

Finally let's do the plot:

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.

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

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)

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:

Int is an alias for the default integer type:

The default floating point type is 64-bit (equivalent of Fortran's double precision or C's double):

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):

But 64-bit floating points are not a subtype of 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).

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

Tuples can be used for multiple variable assignments:

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

Array types, matrices¶

A vector of floats:

A vector of integers:

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:

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):

The size of the array is returned as a tuple:

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

Multiplication and division operators on matrices perform linear algebra operations:

Element-by-element operations are available under .+, .-, .*, ./:

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

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:

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

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:

1. in the declaration of composite types (see above)
2. for qualifying function arguments (see above)
3. to assert that something is of a given type
4. 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:

Here is an example of typed local variables:

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:

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

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):

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:

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:

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.

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

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.

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.

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

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.

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)$:

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.

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$):

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.

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.

Splines¶

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

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).

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.

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).

Verify that VB(0.01) is approximatively -1.01631 and that VG(-0.01, 0.02) is approximatively -0.9959.

Finally, write the function V that is the maximum of VB and VG:

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:

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

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.

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$.

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$).

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$).

Main loop¶

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

Here is a simple example of minimization:

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$).

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.

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
• 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

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