Coursework 2025 - Consumption Saving

Part I : Approximations

Interpolation

Consider the function \(f(x)=\frac{sin(|x|)}{|x|+0.0001}\) defined on \([a,b]\) with \(a=-2\) and \(b=2\).

Consider the grid points \(x_i=a+\frac{i}{N-1}(b-a)\) for \(i\in[0,N]\).

  1. Define in julia a, b, N and the function f . Compute the vector of grid points s (a vector of size N) then the values y (a vector of size N)that f takes at these grid points.
# your code here
  1. Define a function interp(s,y,x) which approximates the original function f at any \(x\in[-2,2]\) but is defined using only s and y (or if you prefer a,b, N and y). (hint: you can use the library Interpolations.jl)
# your code here
  1. On the same plot, represent (s,y) as a scatter plot, the function f and the interpolated values using interp. For the \(x\)-axis choose \([a-0.1, b+0.1]\). How does interp behave outside of \([a,b]\) (i.e. when it extrtapolates). Can you find a way to make it exptrapolate linearly?
# your code here

Discretizing a shock

Consider a random shock \(\epsilon\) that follows a normal law with standard deviation \(\sigma\).

The goal here is just to compute numerically the integral

\[\mathbb{E} f(\epsilon)\]

  1. Define a variable \(\sigma\). Choose a function f so that you can compute the above integral in closed form.
# your code here

Monte-Carlo

Take an integer \(K>0\).

  1. Try to evaluate the integral above by running \(K\) random draws i.e. by using the formula \(f^{MC, K} = \frac{1}{K}\sum_{i=1}^K f(\epsilon_i)\) where \(\epsilon_i\) is a random draw. Ths is the Monte-Carlo approach.
# your code here
  1. Evaluate the standard deviation in \(f^{MC, K}\) by evaluating this function repeatedly. How does the standard deviation depend on \(K\) ?
# your code here

Gauss-Hermite

Discretizing the shock consists in finding real numbers \((x_1, ... x_K)\) and positive real numbers \((w_1, ... w_K)\) such that for a “suitable” function \(f(\epsilon)\) we can write:

\[\mathbb{E} f(\epsilon) \approx \sum_{i=1}^K w_i f(e_i)\]

The Gauss-Hermite quadrature is such a discretization scheme. It is widely documented and there are good Julia libraries for it.

  1. Approximate the expectation above using Gauss-Hermite nodes and weights. What is the error? How does it depend on \(K\) for the function you have chosen?
# your code here

Part II : The Ayiagari Model

In this exercise, we try to solve the Ayiagari model.

The Ayiagari model features a continuum of households indexed by \(i\in[0,1]\), who take wage rate \(\color{red} w\) and interest rate \(\color{red} r\) as given.

Each households \(i\) is hit by a random productivity shock (\(\color{green}e^i_t\)), normally distributed, with variance \(\sigma\) and mean \(1\) so that, at the aggregate, the average productivity is \(1\).

Each household supplies \(\color{green}e^i_t\) units of work, and can decide to consume \(c^i_t\) or save her available income \(a^i_{t}\).

Available income \(a^i_t\) thus follows the law of motion:

\[a^i_{t+1} = {\color{red}w } {\color{green}e^i_{t+1}} + (a^i_{t} - c^i_{t}) {\color{red}r}\]

There is no borrowing so that \(0<c^i_t\leq a^i_{t}\). The household choooses consumption optimally in order to maximize:

\[\mathbb{E}_0 \sum_t \beta^t U(c^i_t)\]

with \(U(x)=\frac{x^{1-\gamma}}{1-\gamma}\) and \(\beta=0.96\), \(\gamma=4.0\).

Aggregate savings are turned into productive capital:

\[K = \int_i a^i_t \]

Note that in equilibrium the distribution of savings is invariant, hence the absence of a time-subscript for aggregate capital.

Also, the ergodicity property implies that the invariant distribution across all agents is identical to the distribution of asset levels reached across time by any single agent \(i\) so that we also have:

\[K = \mathbb{E} [ a^i_t ]\]

Capital is rented by competitive firms that all produce with the same Cobb-Douglas technology. Total production is:

\[Y = K^{\alpha} L^{1-\alpha}\]

with \(L=\overline{L}=1\) and \(\alpha=0.3\).

Wage rate and interest rate are then determined by marginal conditions:

\[{\color{red}r}=\alpha \frac{Y}{K}\] \[{\color{red}w}=(1-\alpha) \frac{Y}{L}\]

From the structure of the whole model, it is clear that capital \(\overline{K}\) determines the whole equilibrium.

Note

Because there is no depreciation, the rental rate determined by the marginal conditions, equals the returns on savings when they are turned into productive capital.

Aggregate production

  1. Create a namedtuple (or a structure) m to hold all the model informations.
# your code here
  1. Create a function rates(K,m)::Tuple{Float64,Float64} that takes aggregate capital and return the rental rate of capital and the wage rate.
# your code here
  1. Choose an initial level of capital \(K_0\) and corresponding rates \(p_0=(r_0,w_0)\) so that \(\beta < \beta r_0 < 1\). Define the corresponding varaibles K_0 and p_0.
Note

The second part of this inequality ensures that wealth distribution is well defined – agents don’t want to accumulate assets until they have unbounded wealth.

Consumption-Savings Problem

We now solve the individual problem given the vector of rates p. This can be done using any solution method (any version of time-iteration for instance). Here we choose to use the classic value-function iteration instead.

Given the formulation of the household problem, we can look for a decision rule for consumption \(c(a)\) with \(c(a)\in]0,a]\). To implement value function iteration, we will also define \(V(a)\).

Note

Since all agents are identical ex ante, we need to solve only one optimization problem. Hence we drop the \(i\) subscripts.

In theory, the state-space is \(]0,\infty]\) but in practice we approximate the system on \(]0,\overline{A}]\) and adjust \(\overline{A}\) until the solution is well defined.

Note

If \(\color{green}e^i_t\) has a compact support, theoretical considerations show that one can choose such an \(\overline{A}\).

  1. Write the Bellman equation
your text there

Approximate space

To represent functions (\(c()\) and \(V()\)) using a finite number of parameters, we discretize the state-space into a finite grid g, made of \(N\) linearly spaced points between \(\epsilon>0\) and \(\overline{A}\) (the role of \(\epsilon>0\) is to avoid undefined behaviour for \(c=0\)).

We use gauss-hermite nodes and weights weights to approximate the shock \(\epsilon\).

  1. Create a namedtuple dis to represent the discretized model (i.e. the grid and the quadrature).
# your code here

Value Function Iteration

  1. Take an initial guess vvec (preferably concave) for the values on the grid and define the function value_update(a, c, vvec, m, p, dis)::Float64 which compute the expression inside the max in the Bellman equation for any state a, any acceptable consumption level c and a continuation value, obtained by interpolating the values vvec at any point outside the grid.
# your code here
  1. Use your preferred method for constrained optimization to compute the value \(c\in]0,a]\) which maximizes the function value_update for a given a of your choice.
Note

This is the crucial step. You can test whether it works for different values of a and debug what’s going by making a plot.

# your code here
  1. Write a function bellman_step(vvec, m, p, dis)::Tuple{Vector, Vector} which takes in a vector representing the value function tomorrow and returns a new value vector and policy vector resulting from the maximization.
# your code here
  1. Write a vfi(m, p, dis) function, which solves the consumption savings problem by value function iteration. It should return the vector of values and the consumption vector.
# your code here
  1. Check that the solution makes sense: plot the solution, the boundaries, check the extrapolation behaviour…
# your code here

Policy Function Iteration

  1. Solve the consumption savings problem using policy function iteration (i.e. using howard improvement steps). Compare execution time with value function iteration.
Note

You might need to define some intermediary functions, like value_step which updates the value at all grid points, for a given vector of policy choices and policy_eval which returns the value vector obtained by iterating value_step many times.

# your code here

Computing the Stable Distribution

Now that we have a consumption rule represented by a vector of consumption values cvec, we would like to compute the corresponding ergodic distribution in order to approximate capital supply \(\int_i a_i\)

With Monte Carlo

  1. Write a function transition(a, cvec, m, p) which returns a new random asset level, from an initial a. The consumption choices are defined using cvec and can be interpolated as before.
# your code here
  1. Propose a way to draw \(L=1000\) random points from the ergodic distribution. Plot the result.
# your code here
  1. Compute the average of these values. It corresponds to the capital supply. What is the standard deviation of this method? (you can evaluate the standard deviation by performing the procedure again)
# your code here

With a Markov Chain

  1. Write a function transition(a, c, m, p, dis)::Vector{Tuple{Float64, Float64}} which returns a vector of asset levels with matching probabilities, obtained from initial level a, with consumption choice c for the various realizations of the discretized shocks (obtained from the quadrature).
# your code here
  1. Use the transition function to define a stochastic matrix P of size \(N\times N\) such that \(P_{ij}\) represents the probability of reaching grid point \(j\) from grid point \(i\). Use the trembling hand method to deal with asset levels that are out of the approximation grid.
Tremblilng hand

When a given asset level \(a\) reached with probability \(w\) falls between two grid points \(a_j\) and \(a_{j+1}\), we consider that it reaches \(a_j\) with probability \(w\frac{a_{j+1}-a}{a_{j+1}-a_j}\) and \(a_{j+1}\) with probability \(w\frac{a_{a-a_j}{a_{j+1}-a_j}\). This trick ensures that the stochastic matrix depends smoothly on the decision rule.

# your code here
  1. Compute the stable distribution \(\mu\) of \(P\). Compute the capital demand.
# your code here

Solve the Model

  1. Using everything you have done so far write a function capital_supply(K, m, p, dis) which computes the rates corresponding to K, solves the consumption problem, computes the ergodic distribution and the corresponding capital supply.
# your code here
  1. Plot capital supply and capital demand (line \(K=K\)) for different levels of capital. Find the level of capital such that market for capital clears.
# your code here