Brock and Mirman (1972)

• The social planner solves the following problem: \[ \begin{gathered} \max _{\left\{c_{t}, k_{t+1}\right\}} \sum_{t=0}^{\infty} \beta^{t} \log \left(c_{t}\right) \\ c_{t}+k_{t+1} \leq A k_{t}^{\alpha} \\ k_{0} \text{ given, } c_{t}, k_{t} \geq 0 \text{ for all } t \end{gathered} \]
• The controls are \(c_{t}\) and \(k_{t+1}\), and the only state is \(k_{t}\). We can reformulate the previous problem as a dynamic programming problem as follows: \[ V(k) = \max_{k'} \left\{ \log \left( A k^{\alpha} - k' \right) + \beta V(k') \right\} \]
• If we solve for capital stock policy rule \(k' = g_k(k)\), we can obtain the policy function for consumption \(c = g_c(k)\) directly from the resource constraint: \[ g_c(k) = A k^{\alpha} - g_k(k) \]

Step 1: Initialization of the parameters

import numpy as np
import matplotlib.pyplot as plt
import time

# Model Parameters
alpha = 0.4
A = 5
beta = 0.9888

# Algorithm Parameters
T = 20  # periods to simulate
tolv = 1e-8  # tolerance

Steady State

k_ss = (1 / (beta * alpha * A)) ** (1 / (alpha - 1))
y_ss = A * k_ss**alpha
c_ss = y_ss - k_ss
ku = (1 / A) ** (1 / (alpha - 1))  # if c = 0

print("k_ss=", k_ss)
print("k_u=", ku)

k_ss= 3.1157606561253433
k_u= 14.62008869106433

Step 2: Creation of a grid for the states

• Create a vector for \(k\), with \(l_{k}\) values i.e. \(k \in [k_1 < k_2 < \ldots < k_{l_k}]\).
• The larger \(l_k\), the better the approximation, but this increases the computing time.
• We can have the grid from \(0\) to \(k_{l_k}\), or around the steady state, e.g. \(\pm 10\%\) of \(k_{ss}\).

lk = 500
k_min = 0.5 * k_ss
k_max = 1.5 * k_ss
k = np.linspace(k_min, k_max, num=lk)

print(k[:5])

[1.55788033 1.56412434 1.57036835 1.57661236 1.58285637]

Step 3: Computation of the return function

• For each value of \(k\) and \(k'\) in the grid we need: \[ \begin{aligned} c &= A k^{\alpha} - k' \\ U(c) &= \log(c) \end{aligned} \]
• \(U\) is of dimension \(g_k \times 1\), where \(g_k = l_k \times l_k\).
• If \(l_k = 2\) , \(U\) is of dimension \(4 \times 1\): \[ \begin{aligned} & \leftarrow k' \longrightarrow \\ & \begin{aligned} &\uparrow \\ & k \\ &\downarrow \end{aligned} \left[ \begin{array}{l} c_{11} \\ c_{12} \\ c_{21} \\ c_{22} \end{array} \right] \\ c_{ij} = A k_{i}^{\alpha} - k_{j} \quad \text{for} \ i=1:l_k,\, j=1:l_k \end{aligned} \]

Step 4: Initialization and computation of the optimal value function

• Our candidate value function is an \(l_k \times 1\) vector: \[ V(k) = \begin{bmatrix} V(k_1) \\ V(k_2) \\ \vdots \\ V(k_{l_k}) \end{bmatrix} \]
• We can initialize \(V_0\) with a vector of zeros of dimension \(l_k\).
• We also need to initialize the contraction mapping operator \(V_1 = T(V_0)\).
• In iteration \(n\), it returns: \[ V_n(k_i) = \max_{k'} \left\{ \begin{bmatrix} U(A k_i - k_1') \\ \vdots \\ U(A k_i - k_{l_k}') \end{bmatrix} + \beta \begin{bmatrix} V_{n-1}(k_1') \\ \vdots \\ V_{n-1}(k_{l_k}') \end{bmatrix} \right\} \]

Value Function Iteration

In each iteration:

• Set \(V_0\) to be the \(V_1\) computed from the previous iteration.
• For each \(k_i\), compute \(U(c) + \beta V_0(k')\) and choose the maximum over \(k'\) and call it \(V_1\) (Use the max operator).
• Compare \(V_0\) with \(V_1\) and, if the difference is sufficiently small, stop (Use the norm of \(V_1 - V_0\)).

# Initialize value function and operator and iterate
V0 = np.ones((lk))
V1 = np.zeros((lk))
start_time = time.time()
while abs(np.linalg.norm(V1 - V0)) > tolv:
    V0 = np.copy(V1)  # Use deepcopy here to avoid shallow copy.
    for i in range(lk):
        vtemp = u[i * lk : (i + 1) * lk] + beta * V0
        V1[i] = np.max(vtemp)

print(V1[:5])
print(f"Time taken: {time.time() - start_time:.2f} seconds")

[138.88421524 138.88686227 138.88949887 138.89212395 138.89474037]
Time taken: 1.57 seconds

Convergence Criteria: Norm Choices

When checking convergence of value function iteration, we often monitor

\[ \|V_{n+1}-V_n\|\le \texttt{tol}. \] Different norms imply different (sometimes stronger/weaker) stopping rules.

1) Sup-norm (\(\ell_\infty\)): the standard for contraction mappings

\[ \|x\|_\infty = \max_i |x_i|. \]

• Interpretation: the largest pointwise change in the value function.
• Why common: Bellman operators are contractions in \(\|\cdot\|_\infty\) under standard conditions.

diff = V1 - V0
err_inf = np.max(np.abs(diff))                 # ||V1 - V0||_inf
# or: err_inf = np.linalg.norm(diff, ord=np.inf)
print(err_inf)

4.4568082557816524e-10

2) Euclidean norm (\(\ell_2\)): average-magnitude change

\[ |x|_2 = \left(\sum_i x_i^2\right)^{1/2}. \]

• Interpretation: overall magnitude of the change.
• Caution: small \(\ell_2\) does not rule out a large error at a single grid point.

diff = V1 - V0
err_2 = np.linalg.norm(diff, ord=2)            # ||V1 - V0||_2
# same as: err_2 = np.sqrt(np.sum(diff**2))
print(err_2)

9.96425562765096e-09

3) \(\ell_1\) norm: total absolute change

\[ |x|_1 = \sum_i |x_i|. \]

• Interpretation: aggregate absolute change across all grid points.
• Caution: depends on grid size (more points → larger sum), so tolerances are not comparable across lk.

diff = V1 - V0
err_1 = np.linalg.norm(diff, ord=1)            # ||V1 - V0||_1
# same as: err_1 = np.sum(np.abs(diff))
print(err_1)

2.2280752887127164e-07

4) Relative norms: scale-free stopping rule

Absolute errors depend on the scale of \(V\). A relative rule is \[ \frac{|V_{n+1}-V_n|}{\max(1,|V_n|)} \le \texttt{tol}. \]

diff = V1 - V0
den = max(1.0, np.linalg.norm(V0, ord=np.inf))
rel_err_inf = np.max(np.abs(diff)) / den
print(rel_err_inf)

3.192300058330576e-12

5) Weighted norm (optional): emphasize important regions

Given weights \(w_i>0\),

\[ |x|_{w,\infty} = \max_i \frac{|x_i|}{w_i}. \]

Useful if some grid points matter more (e.g., near the ergodic distribution).

diff = V1 - V0
w = np.ones_like(diff)          # replace with your weights (must be > 0)
err_winf = np.max(np.abs(diff) / w)
print(err_winf)

4.4568082557816524e-10

Step 5: Computation of the optimal policy functions

• To calculate \(g_k(k)\), for each point in the grid, find the index where \(U(c) + \beta V^*(k')\) takes its maximum value when \(V^*\) is used.
• The element optim[i] gives the index of the \(k'\) maximizer when the initial capital is \(k_i\).

optim = np.zeros(lk, dtype="int") # "int" type for index. 
for i in range(lk):
    vtemp = u[i * lk:(i + 1) * lk] + beta * V0
    optim[i] = np.argmax(vtemp)
print(optim[:5])

[129 129 130 130 131]

Testing our approximation

The analytical solution to the model is:

\[ \begin{gathered} v(k) = (1-\beta)^{-1} \Big[ \ln (A(1-\alpha \beta)) + \frac{\alpha \beta}{1-\alpha \beta} \ln (A \alpha \beta) \Big] + \frac{\alpha}{1-\alpha \beta} \ln (k) \\ k' = g_k(k) = \alpha \beta A k^{\alpha} \\ c = g_c(k) = (1-\alpha \beta) A k^{\alpha} \end{gathered} \]

To test our approximation, we compare the differences between the analytical and numerical solution.

# True value functions
a = (1 / (1 - beta)) * (np.log(A * (1 - alpha * beta)) + alpha * beta * np.log(A * alpha * beta) / (1 - alpha * beta))
b = alpha / (1 - alpha * beta)
Vk = a + b * np.log(k)
dist = np.linalg.norm(Vk - V0)
print(f"Distance between true and approx. value function: {dist:.8f}")

Distance between true and approx. value function: 0.00038753

Plot the true and approximated value functions:

plt.figure(figsize=(6, 4))
plt.plot(k, V0)
plt.plot(k, Vk)
plt.title("True and Approximated value function")
plt.show()

Simulating the model

1. Choose an initial value for \(k_0\). To do this, select a number between 1 and \(I_k\), indicating the position in the initial capital in the grid. Given this position, we use the capital policy function to extract the optimal first period capital.

2. We then create a loop that:

   • given the position of \(k_0\), finds the position of the optimal next period capital that corresponds to it.
   • uses the initial capital and the position of the new optimal capital to calculate the values of output, consumption, and the next period capital, using the optimal policy vectors.

IID Shocks

Step 1: Initialization of the parameters and steady state

# Model Parameters
alpha = 0.4
beta = 0.9888

# Algorithm Parameters
T = 20  # periods to simulate
tolv = 1e-8  # tolerance

Steady State

• To calculate the steady state, we use: \[ \begin{gathered} \frac{1}{c} = \beta \mathbb{E} \frac{1}{c'} (\alpha \theta' {k'}^{\alpha-1}) \\ c + k' = k^{\alpha} \end{gathered} \]
• With \(\mathbb{E}(\theta') = \mathbb{E}(\theta) = 1\) , and \(x = x'\) for all variables, we get: \[ \begin{gathered} k_{ss} = \left(\frac{1}{\beta \alpha}\right)^{\frac{1}{\alpha-1}} \\ c_{ss} = k_{ss}^\alpha - k_{ss} \\ y_{ss} = k_{ss}^\alpha \end{gathered} \]

# steady state
k_ss = (1 / (beta * alpha)) ** (1 / (alpha - 1))
y_ss = k_ss ** alpha
c_ss = y_ss - k_ss

print("k_ss=", k_ss)

k_ss= 0.21311503110303762

Step 2: Creation of a grid for the states

• Create a grid vector for \(\theta\) and \(k\) with \(l_\theta\) and \(l_k\) values: \[ \theta \in [\theta_1, \ldots, \theta_{l_\theta}],\quad k \in [k_1, \ldots, k_{l_k}] \]

lt = 2
theta = np.linspace(1 + 1.625 / 100, 1 - 1.625 / 100, num=lt)
print(theta)

lk = 500
k_min = 0.5 * k_ss
k_max = 1.5 * k_ss
k = np.linspace(k_min, k_max, num=lk)
print(k[:5])

[1.01625 0.98375]
[0.10655752 0.1069846  0.10741168 0.10783877 0.10826585]

Step 3: Computation of the return function

• For each value of \(\theta, k\) and \(k'\), calculate: \[ \begin{gathered} c = \theta k^{\alpha} - k' \\ U(c) = \log(c) \end{gathered} \]

• Create a matrix of dimensions \(g_k \times l_\theta\), where \(g_k = l_k \times l_k\).

gk = lk * lk
c = np.zeros((gk, lt))
for t in range(lt):
    for i in range(lk):
        for j in range(lk):
            c[i * lk + j][t] = theta[t] * k[i] ** alpha - k[j]
            if c[i * lk + j][t] < 0:
                c[i * lk + j][t] = tolv
print(c[:5])

[[0.3084292  0.29515779]
 [0.30800211 0.29473071]
 [0.30757503 0.29430362]
 [0.30714795 0.29387654]
 [0.30672086 0.29344945]]

To calculate the return function, simply take the log of the consumption matrix:

u = np.log(c)
print(u[:5])

[[-1.17626296 -1.22024518]
 [-1.17764863 -1.2216932 ]
 [-1.17903622 -1.22314332]
 [-1.18042574 -1.22459554]
 [-1.18181719 -1.22604987]]

Step 4: Initialization and computation of the optimal value function

• Our candidate solution for the value function is an \(l_k \times l_\theta\) matrix: \[ V(k, \theta) = \begin{bmatrix} V(k_1, \theta_1) & V(k_1, \theta_2) \\ V(k_2, \theta_1) & V(k_2, \theta_2) \\ \vdots & \vdots \\ V(k_{l_k}, \theta_1) & V(k_{l_k}, \theta_2) \end{bmatrix} \]
• Initialize \(V_0\) as zeros of dimension \(l_k \times l_\theta\) and \(V_1 = T(V_0)\). At iteration \(n\), calculate: \[ V_n(k_i, \theta_t) = \max_{k'} \left\{ \begin{bmatrix} U(\theta_t k_i^\alpha - k_1') \\ \vdots \\ U(\theta_t k_i^\alpha - k_{l_k}') \end{bmatrix} + \beta \begin{bmatrix} \sum_{m=1}^2 \frac{1}{2} V_{n-1}(k_1', \theta_m') \\ \vdots \\ \sum_{m=1}^2 \frac{1}{2} V_{n-1}(k_{l_k}', \theta_m') \end{bmatrix} \right\} \]

Value Function Iteration

# Initialize value function and operator and iterate
V0 = np.ones((lk, lt))
V1 = np.zeros((lk, lt))
p = np.array([0.5, 0.5])
start_time = time.time()
while abs(np.linalg.norm(V1 - V0)) > tolv:
    V0 = np.copy(V1)
    for j in range(lt):
        for i in range(lk):
            vtemp = u[i * lk : (i + 1) * lk, j] + beta * V0 @ p
            V1[i, j] = np.max(vtemp)
V0 = np.copy(V1)
print(V0[:5])
print(f"Time taken: {time.time() - start_time:.2f} seconds")

[[-100.6080923  -100.66186288]
 [-100.60544596 -100.65921516]
 [-100.60280894 -100.65657915]
 [-100.60018353 -100.65395331]
 [-100.59756785 -100.65133759]]
Time taken: 5.05 seconds

Step 5: Computation of the optimal policy functions

For each point in the grid, we find the index of the capital maximizer using the max operator.

optim = np.zeros((lk, lt), dtype="int") # "int" type for index. 
for j in range(lt):
    for i in range(lk):
        vtemp = u[i * lk:(i + 1) * lk, j] + beta * V0 @ p
        optim[i, j] = np.argmax(vtemp)
print(optim[:5])

[[135 123]
 [135 123]
 [136 124]
 [137 124]
 [137 125]]

Using this matrix, the policy functions can be calculated as follows:

\[ \begin{gathered} \text{polk} = \begin{bmatrix} g_k(k_1, \theta_1) & g_k(k_1, \theta_2) \\ \vdots & \vdots \\ g_k(k_{l_k}, \theta_1) & g_k(k_{l_k}, \theta_2) \end{bmatrix}, \text{polc} = \begin{bmatrix} g_c(k_1, \theta_1) & g_c(k_1, \theta_2) \\ \vdots & \vdots \\ g_c(k_{l_k}, \theta_1) & g_c(k_{l_k}, \theta_2) \end{bmatrix} \\ \quad \text{and} \quad \text{polc} = \begin{bmatrix} \theta_1 k_1^\alpha & \theta_2 k_1^\alpha \\ \vdots & \vdots \\ \theta_1 k_{l_k}^\alpha & \theta_2 k_{l_k}^\alpha \end{bmatrix} - \text{polk} \end{gathered} \]

polk = k[optim]
print(polk[:5])

[[0.16421389 0.15908888]
 [0.16421389 0.15908888]
 [0.16464097 0.15951596]
 [0.16506806 0.15951596]
 [0.16506806 0.15994304]]

Simulating the model

1. Choose an initial value for \(k_0\) and \(\theta_0\) and use the policy matrix to extract the optimal first-period capital.

2. We then create a loop that:

   • given the initial capital and shock, finds the position of the optimal next period capital.
   • uses the initial capital, the initial shock index, and the next capital to calculate output, consumption, and stores variables for each period. It also generates a new shock index for next period.

Plotting the results

plt.figure(figsize=(6, 6))
plt.plot(range(T), output, label="Output")
plt.plot(range(T), cons, label="Consumption")
plt.plot(range(T), kopt[:T], label="Capital")
plt.title("Simulation of all variables")
plt.xticks(range(0, T + 1, 4))
plt.xlabel("Period")
plt.legend()
plt.show()

Markov Shocks

Markov-Chain Approximation of AR(1) Productivity Shocks

Let the (finite) state space be \(S= \{\theta_1,\dots,\theta_N\}\).

• Markov property (memorylessness). For any \(i,j\), \[ \Pr(\theta_{t+1}=\theta_j \mid \theta_t=\theta_i,\theta_{t-1},\dots) =\Pr(\theta_{t+1}=\theta_j\mid \theta_t=\theta_i)=P_{ij}. \]
• Transition matrix \(P\in\mathbb R^{N\times N}\): \[ P_{ij}\ge 0,\qquad \sum_{j=1}^N P_{ij}=1\ \ \forall i=1,2,...,N. \]

• Stationary (invariant) distribution. A vector \(\pi\in\Delta^{N-1}\) is stationary if \[ \pi^\top=\pi^\top P,\qquad \pi_i\ge 0,\ \sum_i \pi_i=1. \] If the chain is irreducible and aperiodic (ergodic), then \(\pi\) is unique and for any initial distribution \(\mu_0\), \[ \mu_t^\top=\mu_0^\top P^t \to \pi^\top. \]

• How expectations enter DP. If current shock is \(\theta=\theta_i\), \[ \mathbb E\!\left[V(k',\theta')\mid \theta=\theta_i\right]\ \approx\ \sum_{j=1}^N P_{ij}\,V(k',\theta_j). \]

Approximating AR(1) shocks: Tauchen (1986) discretization

• It is standard to discretize log productivity to ensure positivity.

• Assume \[ x_t=\rho x_{t-1}+\varepsilon_t,\qquad \varepsilon_t\sim N(0,\sigma^2),\qquad \theta_t=\exp(x_t). \]

• (i) Choose a grid of \(N\) states for \(x\)

  • Unconditional standard deviation: \[ \sigma_x=\frac{\sigma}{\sqrt{1-\rho^2}}. \]
  • Pick a truncation width \(m>0\) (often \(m=3\)) and set \[ x_1=-m\sigma_x,\quad x_N=m\sigma_x,\quad x_i=x_1+(i-1)\Delta,\quad \Delta=\frac{x_N-x_1}{N-1}. \]

• (ii) Build the transition matrix
  • Define interval cutoffs \[ b_{j-}=x_j-\frac{\Delta}{2},\qquad b_{j+}=x_j+\frac{\Delta}{2}, \] with \(b_{1-}=-\infty\) and \(b_{N+}=+\infty\). Then \[ P_{ij} =\Pr(x_{t+1}\in[b_{j-},b_{j+}]\mid x_t=x_i) =\Phi\!\left(\frac{b_{j+}-\rho x_i}{\sigma}\right)-\Phi\!\left(\frac{b_{j-}-\rho x_i}{\sigma}\right), \] where \(\Phi(\cdot)\) is the standard normal CDF.

Step 1: Initialization of the parameters and steady state

# Model and Algorithm Parameters
delta = 0.1
alf = 0.33
A = 1
beta = 0.9
gamma = 2
tolv = 1e-7

Construction of a Markov Chain

We use the quantecon library to construct a Markov Chain.

• First, we set up the parameter values for the AR(1) process that drives the aggregate productivity shock.
• We then use Tauchen’s method to discretize this process into a finite Markov chain.

# Step 1: Specify the AR(1) process parameters and discretize with Tauchen's method
import quantecon as qe

rho = 0.95  # persistence of the AR(1) process
sigmae = 0.00712  # standard deviation of the error term
lt = 4  # number of grid states for the Markov chain
n_std = 3  # number of standard deviations to cover in the grid

# Use Tauchen's method to approximate the AR(1) with a Markov chain
x_mc = qe.markov.approximation.tauchen(rho=rho, sigma=sigmae, n_std=n_std, n=lt, mu=0)

• Next, we extract the transition matrix, the grid for the state variable representing productivity shocks, compute the values of productivity in levels (by exponentiating), and obtain the stationary distribution over states.

# Step 2: Extract transition matrix and productivity states, and compute stationary distribution
P = x_mc.P  # transition matrix (to match standard notation)
x_grid = x_mc.state_values  # discretized grid for the state variable
teta = np.exp(x_grid)  # productivity in levels, since shock is in logs
invdist = np.array(x_mc.stationary_distributions)  # stationary distribution
print(f"Stationary distribution: {invdist}")
print(f"Productivity states: {teta}")
print(f"Transition matrix: {P[0,:]}")
print(f"Check (sum of first row): {np.sum(P[0,:])}")

Stationary distribution: [[0.05317954 0.44682046 0.44682046 0.05317954]]
Productivity states: [0.93388054 0.97745576 1.02306421 1.07080077]
Transition matrix: [0.99675735 0.00324265 0.         0.        ]
Check (sum of first row): 1.0

Steady State

Eteta = invdist @ teta
k_ss = (A * Eteta * alf * beta / (1 - beta * (1 - delta))) ** (1 / (1 - alf))
y_ss = Eteta * A * (k_ss**alf)
i_ss = delta * k_ss
c_ss = y_ss - i_ss
print(k_ss)

[1.94925376]

Step 2: Creation of a grid for the states

Create a grid for the shock and for \(k\) with \(l_k\) values, i.e. \(k \in [k_1 < \ldots < k_{l_k}]\).

# Grid for capital
lk = 500
k_min = 0.9 * k_ss
k_max = 1.1 * k_ss
k = np.linspace(k_min, k_max, num=lk)
print(k[:5])

[[1.75432839]
 [1.75510965]
 [1.75589092]
 [1.75667218]
 [1.75745344]]

Step 3: Computation of the return function

For each \(\theta, k\) and \(k'\), calculate: \[ \begin{aligned} c = \theta A k^\alpha + (1-\delta)k - k' \\ U(c) = \frac{c^{1-\gamma}}{1-\gamma} \end{aligned} \]

gk = lk * lk
c = np.zeros((gk, lt))
for t in range(lt):
    for i in range(lk):
        for j in range(lk):
            c[i * lk + j][t] = A * teta[t] * k[i] ** alf + (1 - delta) * k[i] - k[j]
            if c[i * lk + j][t] < 0:
                c[i * lk + j][t] = tolv
print(c[0, :])

[0.9487784  1.00123451 1.05613825 1.11360382]

To calculate the return function, take logs if \(\gamma = 1\), or compute \(\frac{c^{1-\gamma}}{1-\gamma}\) otherwise.

if gamma == 1:
    U = np.log(c)
else:
    U = c ** (1 - gamma) / (1 - gamma)
print(U[0, :])

[-1.0539869  -0.99876701 -0.94684574 -0.89798543]

Step 4: Initialization and computation of the optimal value function

Our candidate solution for the value function is an \(l_k \times l_\theta\) matrix.

# Initialize of the value function
V0 = np.ones((lk, lt))
V1 = np.zeros((lk, lt))
start_time = time.time()
while abs(np.linalg.norm(V1 - V0)) > tolv:
    V0 = np.copy(V1)
    for j in range(lt):
        for i in range(lk):
            vtemp = U[i * lk : (i + 1) * lk, j] + beta * V0 @ P[j, :]
            V1[i, j] = np.max(vtemp)

V0 = np.copy(V1)
print(V0[0, :])
end_time = time.time()
print(f"Time taken: {end_time - start_time} seconds")

[-10.5240125   -9.98068286  -9.45989273  -8.96829951]
Time taken: 2.6193039417266846 seconds

Step 5: Computation of the optimal policy functions

optim = np.zeros((lk, lt), dtype="int")
for j in range(lt):
    for i in range(lk):
        vtemp = U[i * lk : (i + 1) * lk, j] + beta * V0 @ P[j, :]
        optim[i, j] = np.argmax(vtemp)
print(optim[0, :])

[ 0 20 42 64]

• For each pair \((\theta_j, k_i)\), optim[i, j] is the index of the next period capital maximizer.
• If \(l_k = 2\) and \(l_\theta = 4\), the policy functions are: \[ \begin{aligned} & \text{polk} = \begin{bmatrix} g_k(k_1, \theta_1) & g_k(k_1, \theta_2) & g_k(k_1, \theta_3) & g_k(k_1, \theta_4) \\ g_k(k_2, \theta_1) & g_k(k_2, \theta_2) & g_k(k_2, \theta_3) & g_k(k_2, \theta_4) \end{bmatrix} \\ & \text{polc} = \begin{bmatrix} \theta_1 k_1^\alpha & \theta_2 k_1^\alpha & \theta_3 k_1^\alpha & \theta_4 k_1^\alpha \\ \theta_1 k_2^\alpha & \theta_2 k_2^\alpha & \theta_3 k_2^\alpha & \theta_4 k_2^\alpha \end{bmatrix} + (1-\delta) \begin{bmatrix} k_1 & k_1 & k_1 & k_1 \\ k_2 & k_2 & k_2 & k_2 \end{bmatrix} - \text{polk} \end{aligned} \]

Simulating the model

1. Choose an initial value for \(k_0\) and \(\theta_0\) and use the policy matrix to extract the optimal first period capital.
2. In a loop:
   • given the initial capital and shock, find the optimal next period capital.
   • use the current values to compute and store output, consumption, and capital, generating a new shock index for next period.

Irreversible Investment

The social planner solves the following problem:

\[ \begin{gathered} \max_{\{c_t, k_{t+1}\}} E_0 \sum_{t=0}^\infty \beta^t \frac{c_t^{1 - \gamma}}{1 - \gamma} \\ c_t + i_t \leq \theta_t A k_t^{\alpha} \\ k_{t+1} = i_t + (1 - \delta) k_t \\ i_t \geq 0 \end{gathered} \]

\(k_0, \theta_0\) given, \(c_t, k_t \geq 0\) for all \(t\).
The code is basically the same as above, except for a small modification in step 2.