Homework 3

Pedro Augusto Januzzi Guerra

My solution to this homework closely followed Ellen's lecture notes and Sergio Ocampo's notes (available in the Holy Grail).

Even though it was not asked, let's first define a Tax Distorted Competitive Equilibrium (TDCE) so that it's easier to solve for the steady state of the problem. I assume the existence of a representative household and firm. For notation simplicity, I do not denote denote variables as being dependend on a specific state, but they are (e.g., $c_t$ should actually be denoted by $c_t(s^t)$). Moreover, I consider the stochastic process for $g_t$ to happen in levels, instead of in log, since this makes the algebra easier.

Definition: A TDCE, given fiscal policy $G = \{\tau_{h_t},\tau_{x_t},g_t\}_{t=0}^\infty$ is a sequence of household's allocations $A^h = \{c_t,x_t,k_{t+1},h_t,l_t\}_{t=0}^\infty$, firm's allocations $A^f = \{K^f_t,L^f_t\}_{t=0}^\infty$ and prices $\{r_t,w_t\}_{t=0}^\infty$ such that:

• Given fiscal policy and prices, household chooses $A^h$ to solve the following maximization problem: \begin{align*} \max_{\hat{A}^h} \hspace{2mm}&\mathbb{E}\sum_{t=0}^\infty \beta^t\{\log(\hat{c}_t) + \psi\log(\hat{l}_t)\}N_t\\[2pt] \text{s.t.}\hspace{2mm} & \hat{c}_t + (1+\tau_{x_t})\hat{x}_t = r_t\hat{k}_t + (1-\tau_{h_t})w_t\hat{h}_t + \kappa_t\\ & N_{t+1}\hat{k}_{t+1} = [(1-\delta)\hat{k}_t+\hat{x}_t]N_t\\ & \hat{h}_t + \hat{l}_t = 1\\ & N_t = (1+\gamma)^t\\ & S_t = PS_{t-1} + Q\epsilon_t \hspace{3mm} \forall S_t \in \{\log(z_t),\tau_{h_t},\tau_{x_t},g_t\}\\ & \hat{c}_t,\hat{x}_t,\hat{h}_t,\hat{l}_t \geq 0 \hspace{2mm}\text{in all states and periods} \end{align*}

• Firms choose $A^f$ to solve: \begin{align*} \max{\hat{A}^f} &\hat{K}^{f^\theta}(z_t\hat{L}^f_t)^{1-\theta}-w_t\hat{L}^f_t - r_t\hat{K}^f_t\\[3pt] \text{s.t.}\hspace{2mm} & \hat{K}^f_t,\hat{L}^f_t \geq 0 \hspace{2mm}\text{in all states and periods} \end{align*}

• All markets clear in all states and periods: \begin{align*} k_t &= K^f_t\\ h_t &= L^f_t\\ N_t(c_t+x_t+g_t) &= K_t^\theta(z_tL_t)^{1-\theta} \end{align*}

• Government budget is balanced node by node: \begin{align*} \tau_{x_t}X_t + \tau_{h_t}w_tH_t = g_t + \kappa_t \end{align*}

in which $X_t$ and $H_t$ denote the aggregate investment and hours worked, respectively.

Now let's calculate the steady state of the model. Once again, I omit the state-dependence for notational simplicity. Let $\tilde{\beta} = \beta*(1+\gamma)$. The FOCs are given by:

\begin{align*} [c_t]:\hspace{2mm}& \frac{\tilde{\beta}^t}{c_t} = \lambda_t\\[2pt] [h_t]:\hspace{2mm}& \frac{\tilde{\beta}^t\psi}{1-h_t} = \lambda_t(1-\tau_{h_t})w_t \\[2pt] [k_{t+1}]:\hspace{2mm}& \lambda_t(1+\tau_{x_t})(1+\gamma)) = \mathbb{E}[\lambda_{t+1}[(1+\tau_{x_{t+1}})(1-\delta)+r_{t+1}]] \end{align*}

Question 2 tells us to assume $\epsilon \sim N(0,\Sigma)$. Thus, $\mathbb{\epsilon} = 0$. In SS, this means that:

$$ \begin{align*} \log(z_1) &= \rho_z \log(z_0) \\ \log(z_2) &= \rho_z \log(z_1) \\ &= \rho_z^2 \log(z_0) \\ &\vdots \\ \log(z_t) &= \rho_z^t\log(z_0) \end{align*} $$

Since $\rho_z \in(0,1)$, it follows that $$\lim_{t\rightarrow \infty} \log(z_t) = 0$$ which in turn implies $z_{ss} = 1$.

Analogously, $g_{ss} = 0$ and $\tau_{x_{ss}} = \tau_{h_{ss}} = 0$.

Now let's compute the Steady State. From $[k_{t+1}]$ in SS, it follows that:

\begin{align*} \frac{\lambda_t}{\lambda_{t+1}} = \frac{(1+\tau_{x_{t+1}})(1-\delta)+r_{t+1}}{(1+\tau_{x_t})(1+\gamma)} \end{align*}

Dividing $[c_t]$ by $[c_{t+1}]$ and using the equation above, we have: \begin{align*} \frac{1}{\tilde{\beta}}\frac{c_{t+1}}{c_t} = \frac{(1+\tau_{x_{t+1}})(1-\delta)+r_{t+1}}{(1+\tau_{x_t})(1+\gamma)} \end{align*}

Dividing $[c_t]$ by $[h_{t}]$, we get the following intratemporal condition: \begin{align*} \frac{1-h_t}{\psi c_t} = \frac{1}{(1-\tau_{h_t})w_t} \end{align*}

The FOCs for the firm's optimization problem + market clearing conditions yield:

\begin{align*} [K^f_t]:\hspace{2mm}& r_t = \theta (z_th_t)^{1-\theta} k_t^{\theta-1}\\ [L^f_t]:\hspace{2mm}& w_t = (1-\theta)z_t^{1-\theta}k_t^\theta h_t^{-\theta} \end{align*}

Using the expression for $r_{t+1}$ that follows from above in the Euler Equation yields:

\begin{align*} \frac{1}{\tilde{\beta}}\frac{c_{t+1}}{c_t} = \frac{(1+\tau_{x_{t+1}})(1-\delta)+\theta (z_{t+1}h_{t+1})^{1-\theta} k_{t+1}^{\theta-1}}{(1+\tau_{x_t})(1+\gamma)} \end{align*}

In SS, the expression above becomes: \begin{align*} \frac{1}{\tilde{\beta}} &= \frac{(1-\delta)+\theta (h_{ss})^{1-\theta} k_{ss}^{\theta-1}}{(1+\gamma)}\\ \Longrightarrow \frac{1}{\beta} &= (1-\delta)+\theta (h_{ss})^{1-\theta} k_{ss}^{\theta-1}\\ \Longrightarrow \frac{1}{\theta}\Bigg(\frac{1}{\beta} - 1 + \delta\Bigg) &= (h_{ss})^{1-\theta} k_{ss}^{\theta-1}\\ \Longrightarrow k_{ss}^{1-\theta}\frac{1}{\theta}\Bigg[\frac{1}{\beta} - 1 + \delta\Bigg] &= h_{ss}^{1-\theta} \end{align*}

Using the expression for $w_{t+1}$ that follows from above in the intratemporal condition yields:

\begin{align*} \frac{1-h_t}{\psi c_t} = \frac{1}{(1-\tau_{h_t})(1-\theta)z_t^{1-\theta}k_t^\theta h_t^{-\theta}} \end{align*}

In SS, the expression above becomes: \begin{align*} \frac{1-h_{ss}}{\psi c_{ss}} = \frac{1}{(1-\theta)k_{ss}^\theta h_{ss}^{-\theta}} \end{align*}

Note that in SS, market clearing for goods imply: \begin{align*} c_{ss} + x_{ss} = k_{ss}^\theta h_{ss}^{1-\theta} \end{align*}

But from the Law of Motion in SS: \begin{align*} (\delta+\gamma) k_{ss} = x_{ss} \end{align*}

Then, $$c_{ss} + (\delta+\gamma) k_{ss} = k_{ss}^\theta h_{ss}^{1-\theta}$$

Let's summarize what we have so far:

\begin{align*} \text{Euler Equation in SS:}\hspace{2mm} & k_{ss}\Bigg(\frac{1}{\theta}\Big[\frac{1}{\beta} - 1 + \delta\Big]\Bigg)^{\frac{1}{1-\theta}} = h_{ss}\\[8pt] \text{Intratemporal condition in SS:}\hspace{2mm} & \frac{1-h_{ss}}{\psi c_{ss}} = \frac{1}{(1-\theta)k_{ss}^\theta h_{ss}^{-\theta}}\\[10pt] \text{Market clearing for goods in SS:}\hspace{2mm} & c_{ss} + (\delta+\gamma) k_{ss} = k_{ss}^\theta h_{ss}^{1-\theta} \end{align*}

This is a system of 3 equations and 3 unkowns, so we can solve for each of the unkowns.

Isolate consumption in the market clearing to get $$c_{ss} = k_{ss}^\theta h_{ss}^{1-\theta} - (\delta+\gamma) k_{ss}$$

Let $\Lambda = \Bigg(\frac{1}{\theta}\Big[\frac{1}{\beta} - 1 + \delta\Big]\Bigg)^{\frac{1}{1-\theta}}$. Substitute the expression for $h_{ss}$ and $c_{ss} in the intratemporal condition and solve for $k_{ss}$ to get: $$k_{ss} = \frac{1-\theta}{\psi(\Lambda - (\delta+\gamma)\Lambda^\theta)+(1-\theta)\Lambda}$$

It clearly follows that for the other conditions that:

\begin{align*} h_{ss} &= k_{ss}\Lambda\\[2pt] c_{ss} &= k_{ss}^\theta h_{ss}^{1-\theta} - (\delta+\gamma) k_{ss} \end{align*}

---

Question 1

# Packages used
using ForwardDiff, LinearAlgebra, Plots, Distributions

# Parameters (chosen by me)
β = 0.95; # discount factor
δ = 0.05; # depreciation rate
θ = 0.33; # output elasticity of capital
γ = 0.01; # population growth
ψ = 1e-2; # relative preference for leisure over consumption

pers = [0.95, 0.99 , 0.99, 0.99]; # persistence for the the stochastic processes: log z, taux, tauh, g

# Steady state values for stochastic processes:
z_ss = 0; # log(z_ss) = 0
τ_xss = 0;
τ_hss = 0;
g_ss = 0;

Create a function to compute the SS:

function ss(θ,β,δ,γ,ψ)

    Λ = ((1/θ)*((1/β)-1+δ))^(1/(1-θ)); # parameter Lambda
    
    k = (1-θ)/((ψ*(Λ - (δ+γ)*Λ^θ))+(1-θ)*Λ); # ss level of capital
    h = k*Λ; # ss level of hours worked
    l = 1 - h; # ss level of leisure
    c = k^θ*h^(1-θ) - (δ+γ)*k; # ss level of consumption
    x = k*(γ+δ); # ss level of investment

    return k,h,l,c,x
end

ss (generic function with 1 method)

Compute SS variables:

k_ss, h_ss, l_ss, c_ss, x_ss = ss(θ,β,δ,γ,ψ);

Create function to compute the utility:

function util(y)
    
    # y is a vector containing the following variables:
    k = y[1]; # individual current capital holdings (individual state)
    z = exp(y[2]); # current shock (aggregate state)
    τ_x = y[3]; # investment taxation (aggregate state)
    τ_h = y[4]; # labor income taxation (aggregate state)
    g = y[5]; # government expenditures (aggregate state)
    K = y[6]; # aggregate capital holdings, which is equal to k, but the HH doesn't know (aggregate state)
    X = y[7]; # law of motion for aggregate X (aggregate state)
    # agg investment is the second aggregate state, since it matters for transfers
    H = y[8]; # aggregate labor, which is equal to h, but the HH doesn't know (aggregate state)
    x = y[9]; # individual investment (control variable)
    # actually, I'll use x as the control variable, but by choosing x I'm already choosing k' and vice-versa.
    h = y[10]; # individual labor choice (control variable)

    # Compute prices using the FOCs of firm's maximization problem:
    w = (1-θ)*(z^(1-θ))*(K^θ)*(H^(-θ)); # wages
    r = θ*((z*H)^(1-θ))*(K^(θ-1)); # interest rate

    # Compute aggregate Kprime
    Kprime = X/(1+γ) + ((1-δ)/(1+γ)*K);

    # Compute individual kprime
    kprime = x/(1+γ) + ((1-δ)/(1+γ)*k);

    # Compute transfers using the government budget balancedness condition:
    kappa = τ_x*X + τ_h*w*H - g;

    # Compute consumption using the budget constraint:
    c = r*k + (1-τ_h)*w*h + kappa - (1+τ_x)*x;

    utility = log(c) + ψ*log(1-h);

    return utility
end

util (generic function with 3 methods)

Compute the utility level in the SS:

#= Define a vector with all the equilibrium variables. Note that the aggregate variables are equal to individual 
(but the agent doesn't know: big K, little k trick). Note also that in SS, prime variables are equal to current variables.=#
y = [k_ss,z_ss,τ_xss,τ_hss,g_ss,k_ss,x_ss,h_ss,x_ss,h_ss];

# Compute utility in SS:
util_ss = util(y);

Regarding the adaptation to LQ approximation, what should we know?

Well, first notice that the initial goal is the same as before: approximate our return function by a quadratic one and our system of constraints by a linear system. Using the same notation as in Homework 2, the approximation for the return function will be:

$$ r(X,u) \approx \begin{bmatrix} X^T & u^T \end{bmatrix} \begin{bmatrix} Q_{11} & Q_{12}^T \\ Q_{12} & Q_{22} \end{bmatrix} \begin{bmatrix} X \\ u \end{bmatrix} = X^TQX + u^TRu + 2X^TWu $$

where $Q_{11} = \bar{r}-\bar{T}^T\bar{J} + \frac{1}{2}\bar{T}^T\bar{H}\bar{T}$, $\hspace{3mm}$ $Q_{12}=\frac{1}{2}(\bar{J}-\bar{H}\bar{T})$, $\hspace{3mm}$ and $Q_{22} = \frac{1}{2}\bar{H}$.

The main difference relative to HW 2 is that, since now we don't have a planner's problem, the structure and dimensions of the matrices characterizing the solution to the problem will be different. Hence, in order to apply the same methods as in HW 2, we will have to adapt our notation.

First, note that now the vector $X$ is given by $\begin{bmatrix} X_1 & X_2 & X_3\end{bmatrix}'$, where $X_1$ are the individual states (and includes the constant), $X_2$ are the aggregate states with known laws of motion, while $X_3$ are the aggregate states with laws of motion that are unknown, but need to be computed in equilibrium. In our model, the only individual state is $k$, so $X_1 = \begin{bmatrix} 1 & k_t\end{bmatrix}'$; the aggregate states with known laws of motion are $\begin{bmatrix} \log z_t & \tau_{h_t} &\tau_{x_t} & g_t\end{bmatrix}'$; and $X_3$ is composed by $K$, $H$, and $X$ (note that instead of $X$, we could write $K'$ since both are related to each other; in fact, in the code I'll work with $K'$). Note that investment matters in this problem since we need them to compute total transfers using the government budget balance condition.

Lastly, note from the paragraph above then that $X$ will be a vector of dimensions $9$ x $1$ and $u$ a vecor of dimensions $2$ x $1$, since we will work with two control variables - $x_t$ and $h_t$. Note that, with those, we can compute all other control varaibles: using the time constraint condition, we compute $l_t$; with the law of motion, we compute $k_{t+1}$; with budget constraint, we compute $c_t$. This means that matrix $Q$ will be $9$ x $9$, $W$ will be $9$ x $2$, and $R$ will be $2$ x $2$.

Create a function to get the Jacobian and Hessian matrices:

function getJacobHess(y,r)

    J = ForwardDiff.gradient(r,y) 
    H = ForwardDiff.hessian(r,y)

    # Round the values in the Jacobian and Hessian to 5 decimal places
    #J = round.(J, digits=5);
    #H = round.(H, digits=5);    

    return J,H
end

getJacobHess (generic function with 1 method)

Create function to compute matrices for LQ approximation:

function matrices(J,H,T,util_ss) # J = Jacobian, H = Hessian, T = vector of variables in SS, util_ss = utility in SS

    Q11 = util_ss - T'*J + 0.5*T'*H*T; # this is Q in my notation above

    Q22 = 0.5*H; # this is R in my notation above

    Q12 = 0.5*(J-H*T);

    # Create a block matrix
    Qaux = [Q11 Q12'; Q12 Q22];

    s = size(Qaux)[1]; # get the number of rows of the block matrix above (since Qaux is symmetrics, s = number of cols)
    # in our case, it will be 11 x 11

    # Partition Q
    Q = Qaux[1:(s-2),1:(s-2)]; # rows 1 to 9, cols 1 to 9: matrix 9x9
    W = Qaux[1:(s-2),(s-1):s]; # rows 1 to 9, cols 10 to 11: matrix 9 x 2
    R = Qaux[(s-1):s,(s-1):s]; # rows and cols 10 to 11: matrix 2 x 2

    return Q,R,W
end

matrices (generic function with 1 method)

Compute the Jacobian and Hessian matrices:

J,H = getJacobHess(y,util);

Compute matrices Q,R,W:

Q,R,W = matrices(J,H,y,util_ss);

Map $Q$ into an undiscounted problem without cross-products in the objective function:

# t stands for tilde
Qt = Q-W*inv(R)*W';

Now, let's move to the matrices that are related to the approximation to the constraints of the problem. Following Ellen's notation, the real system could be written as $$X_{t+1} = g(X_t,u_t,\epsilon_{t+1}) \hspace{3mm}\forall t$$

We will search for matrices $A$, $B$, and $C$ such that: $$g(X_t,u_t,\epsilon_{t+1}) \approx AX_t + Bu_t + C\epsilon_{t+1}$$

That is,

$$ \begin{bmatrix} X_1 \\ X_2 \\ X_3\end{bmatrix}_{t+1} = \begin{bmatrix} A_{11} & A_{12} & A_{13} \\ 0 & A_{21} & A_{23} \\ 0 & A_{32} & A_{33}\end{bmatrix}\begin{bmatrix} X_1 \\ X_2 \\ X_3\end{bmatrix}_{t} + \begin{bmatrix} B_1 \\ 0 \\ 0\end{bmatrix}u_t + C\begin{bmatrix} \epsilon_1 \\ \epsilon_2 \\ \epsilon_3\end{bmatrix}_{t+1} $$

where all matrices above are known, except $A_{32}$ and $A_{33}$, since those are related to the aggregate states with unkown laws of motion.

Following the suggestion in Ellen's notes, "we can stack the states with known laws of motion into a new vector $y_t = \begin{bmatrix}X_{1t} & X_{2t} \end{bmatrix}'$ and work with the upper partition" of the matrix above as follows:

$$ y_{t+1} = A_yy_t + B_yu_t + A_zX_{3t} $$

where $A_y$ is the upper left partition of $A$, $B_y$ is the upper partition of $B$, and $A_z$ the upper right partition of $A$.

Note now that our system $y_{t+1}$ will represent the following equations: \begin{align*} k_{t+1} &= \frac{1-\delta}{1+\gamma}k_t + \frac{1}{1+\gamma}x_t\\ \log z_{t+1} &= \rho_z\log z_t + \varepsilon_{z_{t+1}}\\ \tau_{x_{t+1}} &= \rho_x\tau_{x_{t}} + \varepsilon_{x_{t+1}}\\ \tau_{h_{t+1}} &= \rho_h\tau_{h_{t}} + \varepsilon_{h_{t+1}}\\ g_{t+1} &= \rho_g g_{t} + \varepsilon_{g_{t+1}} \end{align*}

Besides those, we also have the constant $1$ at the beginning of $y_{t+1}$. Therefore, our matrices $A_y$, $A_z$, and $B_y$ are:

\begin{align*} A_y = \begin{bmatrix} 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & \frac{1-\delta}{1+\gamma} & 0 & 0 & 0 & 0 \\ 0 & 0 & \rho_z & 0 & 0 & 0 \\ 0 & 0 & 0 & \rho_x & 0 & 0 \\ 0 & 0 & 0 & 0 & \rho_h & 0 \\ 0 & 0 & 0 & 0 & 0 & \rho_g\end{bmatrix} \end{align*}

\begin{align*} A_z = \begin{bmatrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix} \end{align*}

\begin{align*} B_y = \begin{bmatrix} 0 & 0 \\ \frac{1}{1+\gamma} & 0 \\ 0 & 0 \\ 0 & 0 \\ 0 & 0 \\ 0 & 0 \end{bmatrix} \end{align*}

PS: Note that the full matrix $A$ is given by: \begin{align*} A = \begin{bmatrix} 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & \frac{1-\delta}{1+\gamma} & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & \rho_z & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & \rho_x & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & \rho_h & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & \rho_g & 0 & 0 & 0 \\ 0 & 0 & ? & ? & ? & ? & ? & ? & ?\\0 & 0 & ? & ? & ? & ? & ? & ? & ?\\ 0 & 0 & ? & ? & ? & ? & ? & ? & ?\end{bmatrix} \end{align*}

Computing $A_y$, $A_z$, and $B_y$:

Ay = zeros(6,6);
Ay[diagind(Ay)] = [1, (1-δ)/(1+γ), pers[1], pers[2], pers[3], pers[4]]; # changing the elements in the diagonal

Az = zeros(6,3);

By = zeros(6,2);
By[2,1] = 1/(1+γ); # changing the element in row-2, col-1

Since we are working with partitions of $A$ and $B$, let's get the correspondent partitions of $W$ and $Q$:

Wy = W[1:6,:]; # recall that W is 9x2, so I'm getting now a matrix 6 x 2
Wz = W[7:9,:]; # matrix 3 x 2

Qyt = Qt[1:6,1:6]; # note X1 (including constant) and X2 stacked will be 6x1, so now the correspondent Qyt has to be 6 x 6
Qzt = Qt[1:6,7:9]; # these terms in columns 7 to 9 correspond to those who multiply X3

Let's redefine the matrices to get rid of the cross products, $\gamma$ and $\beta$:

Ayt = sqrt((1+γ)*β)*(Ay - By*inv(R)*Wy'); # t stands for tilde
Azt = sqrt((1+γ)*β)*(Az - By*inv(R)*Wz');
Byt = sqrt((1+γ)*β)*By;

Now we impose market clearing conditions, which are given by:

$$ X_{3t} = \Theta\begin{bmatrix} X_{1t} & X_{2t}\end{bmatrix}' + \Psi u_t $$

In our model, the market clearing conditions for the aggregate states are:

\begin{align*} K_t = k_t\\ X_t = x_t\\ H_t = h_t \end{align*}

So, what we want is to find matrices $\Theta$ and $\Psi$ such that:

\begin{align*} \begin{bmatrix} K \\ X \\ H\end{bmatrix}_t = \Theta\begin{bmatrix} 1\\ k \\ \log z \\ \tau_x \\ \tau_h \\ g\end{bmatrix}_t + \Psi\begin{bmatrix} x \\ h\end{bmatrix}_t \end{align*}

Clealy, in order to satisfy the market clearing conditions written above, we need $\Theta$ and $\Psi$ to satisfy:

\begin{align*} \Theta = \begin{bmatrix} 0 & 1 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 \end{bmatrix} \end{align*}

\begin{align*} \Psi = \begin{bmatrix} 0 & 0 \\ 1 & 0 \\ 0 & 1 \end{bmatrix} \end{align*}

Compute $\Theta$ and $\Psi$:

Θ = zeros(3,6);
Θ[1,2] = 1;
Ψ = [0 0; 1 0; 0 1];

Map the original problem to the undiscounted one:

Θt = inv(I(3) + Ψ*inv(R)*Wz')*(Θ - Ψ*inv(R)*Wy');
Ψt = inv(I(3) + Ψ*inv(R)*Wz')*Ψ;

This yields the following system:

$$ \tilde{X}_{3t} = \tilde{\Theta}\begin{bmatrix} \tilde{X}_{1t} & \tilde{X}_{2t}\end{bmatrix}' + \tilde{\Psi} \tilde{u}_t $$

From the notes:

"Plugging the expression for $\tilde{X}_{3t}$ equation into the first order conditions and doing some simple algebra yields a system just like before:

$$ \begin{bmatrix} y \\ \lambda \end{bmatrix}_t = \begin{bmatrix} \hat{A}^{-1} \hat{Q} \hat{A}^{-1} & \hat{A}^{-1} \hat{B} \hat{R}^{-1} \hat{B}_y' \\ \hat{Q} \hat{A}^{-1} \hat{Q} \hat{B} \hat{R}^{-1} \hat{B}_y' \hat{A} + \hat{A} \end{bmatrix} \begin{bmatrix} y \\ \lambda \end{bmatrix}_{t+1} = H \begin{bmatrix} y \\ \lambda \end{bmatrix}_{t+1} $$

where

$$ \hat{A} = \tilde{A}_y + \tilde{A}_z \hat{\Theta}, \quad \hat{Q} = \tilde{Q}_y + \tilde{Q}_z \hat{\Theta}, \quad \hat{B} = \tilde{B}_y + \tilde{A}_z \hat{\Psi}, \quad \bar{A} = \tilde{A}_y - \tilde{B}_y \hat{R}^{-1} \hat{\Psi}' \tilde{Q}_z $$

The matrix $H$ is analogous to Vaughan's Hamiltonian matrix. Note that if $\hat{\Theta} = 0$ and $\hat{\Psi} = 0$, $H$ is Vaughan's matrix with eigenvalues that come in reciprocal pairs. When there are distortions, we will not find this pattern but the numerics are just as simple."

Compute $\hat{A}$, $\hat{Q}$, $\hat{B}$ and $\bar{A}$:

Ah = Ayt + Azt*Θt; # h stands for hat
Qh = Qyt + Qzt*Θt;
Bh = Byt + Azt*Ψt;
Ab = Ayt - Byt*inv(R)*Ψt'*Qzt'; # b stands for bar

Compute the Hamiltonian matrix:

a11 = inv(Ah);
a12 = inv(Ah)*Bh*inv(R)*Byt';
a21 = Qh*inv(Ah);
a22 = Qh*inv(Ah)*Bh*inv(R)*Byt' + Ab';

ha = [a11 a12; a21 a22];

Get the eigenvalues and eigenvectors, sort the eigenvectors in descending order and get the spectral decomposotion:

# Get the eigenvalues and eigenvectors of the Hamiltonian matrix:
egval,egvec = eigen(ha); 

# Sort them in descending order
egvec = egvec[:,sortperm(egval,rev=true)];

# Spectral decomposition
s = Int(size(egvec)[1]/2);
V_11 = egvec[1:s, 1:s];
V_12 = egvec[1:s, (s+1):(2*s)];
V_21 = egvec[(s+1):(2*s), 1:s];
V_22 = egvec[(s+1):(2*s), (s+1):(2*s)];

Compute $P$:

P = V21*inv(V11); # solution for P

Compute $F$ for the discounted and undiscounted problem:

# Undiscounted:
Ft = inv(R + Byt'*P*Bh)*Byt'*P*Ah;

# Map to the discounted problem:
F = inv(R+Wz'*Ψ)*(R*Ft + Wy' + Wz'*Θ);

Compute the solution for $x$ and $h$:

x_sol,h_sol = -F*[1 k_ss z_ss τ_xss τ_hss g_ss]';

Testing whether this solution is the same as the steady state levels:

isapprox([x_sol, h_sol], [x_ss, h_ss], atol=1e-6)

true

Create a function to compute the value function:

function vf(kgrid,vec,P) #n number of shocks

    z = vec[1];
    τx = vec[2];
    τh = vec[3];
    g = vec[4];

    nk = length(kgrid);
    value_function = zeros(nk);

    for i = 1:nk
        value_function[i] = [1 , k_grid[i] , z_ss , τ_xss , τ_hss , g_ss]' * P * [1 , k_grid[i] , z_ss , τ_xss , τ_hss , g_ss];
    end

    vf = value_function;
    
    return vf
end

vf (generic function with 4 methods)

Create a grid for capital:

# Create the grid for capital
nk = 500; # number of grid points
lb = 0.7*k_ss; # lower bound around 30% below the SS level
ub = 1.3*k_ss; # upper bound around 30% above the SS level
k_grid = LinRange(lb,ub,nk);

Compute the value function:

# compute vector for ss level of stochastic variables
ss_vec = [z_ss τ_xss τ_hss g_ss]; 

# Compute the value function:
value_function = vf(k_grid,ss_vec,P);

Create a function to compute the policy function for investment:

function x_pol(kgrid,vec,F)

    z = vec[1];
    τx = vec[2];
    τh = vec[3];
    g = vec[4];

    nk = length(kgrid);
    x_function = zeros(nk);

    for i = 1:nk
        x_function[i] = (-F*[1 kgrid[i] z τx τh g]')[1];
    end

    return x_function
end

x_pol (generic function with 1 method)

Compute policy function for investment:

x_pf = x_pol(k_grid,ss_vec,F);

Compute policy function for capital:

k_pf = zeros(nk);

k_pf .= x_pf .+ ((1-δ)/(1+γ)).*k_grid;

Plot graph for the value function and policy function for capital:

plot1 = plot(k_grid, value_function, legend=false,title="Value Function", xlabel="Capital Values",linewidth=3);
plot2 = plot(k_grid,k_pf,legend=false,title="Policy Function for Capital",xlabel="Capital Values",linewidth=3);

plot(plot1,plot2,layout=(1,2),size = (1200, 600),bottom_margin = 10Plots.mm)

---

Question 2

First, note that all processes have 0 mean. There is no restriction on the variance, so I'll arbitrarily choose them. To make things interesting, I'll choose different variance parameters for all of them.

Choosing variances:

σ_z = 0.04; # TFP shocks
σ_h = 0.002; # labor income
σ_x = 0.005; # investment
σ_g = 0.01; # government expenditures

variance = [σ_z,σ_x,σ_h,σ_g];

Create a function that receives persistence parameters, variance, steady state levels, number of periods and computes a sequence of observations for our stochastic variables:

function simulation(p,v,ss,t)

    # Get the variances for ϵ
    σz = v[1];
    σx = v[2];
    σh = v[3];
    σg = v[4];

    # Get the persistence parameters
    ρz = p[1];
    ρx = p[2];
    ρh = p[3];
    ρg = p[4];

    # Get the SS variables for the parameters
    zss = ss[1];
    τxss = ss[2];
    τhss = ss[3];
    gss = ss[4];

    # Generate random shocks for all processes (mean zero and individual variance)
    ϵ_zt = rand(Normal(0,sqrt(σz)),t); # epsilon for tfp shocks
    ϵ_xt = rand(Normal(0,sqrt(σx)),t); # epsilon for τ_x
    ϵ_ht = rand(Normal(0,sqrt(σh)),t); # epsilon for τ_h
    ϵ_gt = rand(Normal(0,sqrt(σg)),t); # epsilon for g

    # Create empty vectors to store the sequence of t periods
    z_t = zeros(t);
    τx_t = zeros(t);
    τh_t = zeros(t);
    g_t = zeros(t);

    # Set initial values equal to the steady state
    z_t[1] = zss;
    τx_t[1] = τxss;
    τh_t[1] = τhss;
    g_t[1] = gss;

    # Compute the entire sequence
    for i = 2:t
        z_t[i] = ρz*z_t[i-1] + ϵ_zt[i];
        τx_t[i] = ρx*τx_t[i-1] + ϵ_xt[i];
        τh_t[i] = ρh*τh_t[i-1] + ϵ_ht[i];
        g_t[i] = ρg*g_t[i-1] + ϵ_gt[i];
    end

    return z_t,τx_t,τh_t,g_t # z_t here is actually in logs
end

simulation (generic function with 1 method)

The intuition behind the previous function is to simulate shocks in an economy that was initially at the steady-state.

Drawing observations for 100 periods:

t = 100; # number of periods
sss = [z_ss,τ_xss,τ_hss,g_ss]; # sss = steady state for stochastic (variables)

z_sim, τx_sim, τh_sim, g_sim = simulation(pers,variance,sss,t);

Let's plot the simulations:

p1 = plot(1:t,z_sim,legend=false,xlabel="Periods",title="TFP shocks");
p2 = plot(1:t,τx_sim,legend=false,xlabel="Periods",title="Investment taxes");
p3 = plot(1:t,τh_sim,legend=false,xlabel="Periods",title="Labor income taxes");
p4 = plot(1:t,g_sim,legend=false,xlabel="Periods",title="Government expenditures");

plot(p1,p2,p3,p4,layout=(2,2),size = (1200, 600),bottom_margin = 10Plots.mm)

Now, recall from our HW 2 that (adapting to our current notation): $$u_t = -Fy_t$$

where $u_t =\begin{bmatrix} x_t & h_t\end{bmatrix}$'.

Moreover, recall from question 1 that $y_t = \begin{bmatrix} 1 & k_t & z_t & \tau_{x_t} & \tau_{h_t} & g_t\end{bmatrix}'$.

What we'll do next is to simulate a sequence for $x_t$ and $h_t$ using this rule. Given that we have current capital holdings and investment choice, we can recover the policy function for capital using the law of motion. With those, we can construct similar time series for the other variables.

Construct time series for investment and hours worked:

# Create empty vectors for investment, hours and capital
x_sim = zeros(t);
h_sim = zeros(t);
k_sim = zeros(t);

# Store first entry with steady state values
x_sim[1] = x_ss;
h_sim[1] = h_ss;
k_sim[1] = k_ss;

# PS: here I do not need to take exp(z) because of the way I defined my functions at the beginning of question 1.

for i = 1:t
    if i == 1 
        k_sim[i+1] = ((1-δ)/(1+γ)*k_sim[i]) + (x_sim[i]/(1+γ)); # compute capital for tomorrow using the law of motion
    elseif i > 1 && i < t
        x_sim[i],h_sim[i] = -F*[1 k_sim[i] z_sim[i] τx_sim[i] τh_sim[i] g_sim[i]]'; # compute x and h using the policy function
        k_sim[i+1] = ((1-δ)/(1+γ)*k_sim[i]) + (x_sim[i]/(1+γ));
    elseif i == t
        x_sim[i],h_sim[i] = -F*[1 k_sim[i] z_sim[i] τx_sim[i] τh_sim[i] g_sim[i]]';
    end
end

Let's compute a series for prices using the firm's FOCs:

z_levels = zeros(t);
z_levels = [exp(z_sim[i]) for i=1:t]; # take the exponential to work with the process of z in levels, not in logs

r_sim = zeros(t);
w_sim = zeros(t);

for i = 1:t
    r_sim[i] = θ*((z_levels[i]*h_sim[i])^(1-θ))*(k_sim[i]^(θ-1));
    w_sim[i] = (1-θ)*(z_levels[i]^(1-θ))*(k_sim[i]^θ)*(h_sim[i]^(-θ));
end

Now, let's compute a sequence of transfers using the government budget balance condition.

κ_sim = zeros(t);

κ_sim = [τx_sim[i]*x_sim[i] + τh_sim[i]*w_sim[i]*h_sim[i] - g_sim[i] for i=1:t];

Lastly, let's compute the time series for consumption and leisure:

c_sim = zeros(t);
l_sim = zeros(t);

# Sequence for leisure:
l_sim = [1-h_sim[i] for i = 1:t];

# Sequence for consumption
for i = 1:t
    c_sim[i] = r_sim[i]*k_sim[i] + (1-τh_sim[i])*w_sim[i]*h_sim[i] - (1+τx_sim[i])*x_sim[i] - κ_sim[i];
end

Let's compute the time series for consumption, capital, investment, hours worked, and leisure:

p1 = plot(1:t,c_sim,legend=false,xlabel="Periods",title="Consumption");
p2 = plot(1:t,k_sim,legend=false,xlabel="Periods",title="Capital");
p3 = plot(1:t,x_sim,legend=false,xlabel="Periods",title="Investment");
p4 = plot(1:t,h_sim,legend=false,xlabel="Periods",title="Hours worked");
p5 = plot(1:t,l_sim,legend=false,xlabel="Periods",title="Leisure time");

plot(p1,p2,p3,p4,p5,layout=(3,2),size = (1300, 700),bottom_margin = 10Plots.mm)

Once we have generated the main time series, now we can think about other time series. For instance, the question asks us to construct time series for dividends, accounting profits, and stock valuations. Let's proceed with this analysis.

Dividends are defined as the firm's profits that aren't used for investment.

Compute dividends:

d_sim = zeros(t); # dividends
for i=1:t
    d_sim[i] = (k_sim[i]^θ)*((z_levels[i]*h_sim[i])^(1-θ)) - w_sim[i]*h_sim[i] - (1+τx_sim[i])*x_sim[i];
end

Plot the time series for dividends:

plot(1:t,d_sim,legend=false,xlabel="Periods",title="Dividends")

Accounting profits will be equal to dividends plus capital replacement.

Compute accounting profits (corporate profits):

ap_sim = zeros(t-1); # accounting profits
ap_sim = [d_sim[i]+k_sim[i+1]-k_sim[i] for i=1:(t-1)];

Plot accounting profits:

plot(1:(t-1),ap_sim,legend=false,xlabel="Periods",title="Accounting Profits")

We get stock valuations from Tobin's Q. In the model, $Q = 1+\tau_{x_t}$. From Tobin's Q, we have:

\begin{align*} Q_t &= \frac{v_t}{k_t}\\ \Longrightarrow v_t &= (1+\tau_{x_t})k_t \end{align*}

Compute stock valuations:

sv_sim = zeros(t); # stock valuations
sv_sim = [(1+τx_sim[i])*k_sim[i] for i=1:t];

Plot stock valuations:

plot(1:t,sv_sim,legend=false,xlabel="Periods",title="Stock Valuations")

Another interesting statistics is GDP composition. This is useful to understand what is the contribution of each component of GDP during the business cycle. We know that output $Y_t$ in our environment is given by: $$Y_t = N_t(c_t+x_t+g_t)$$

Compute output:

y_sim = zeros(t);

y_sim = [((1+γ)^i)*(c_sim[i] + x_sim[i] + g_sim[i]) for i=1:t];

Plot GDP with each of its components throughout time:

plot(1:t,y_sim,label="GDP",xlabel="Periods",title="GDP Decomposition")
plot!(1:t,c_sim,label="Consumption")
plot!(1:t,x_sim,label="Investment")
plot!(1:t,g_sim,label="Government Expenditures")

Compute the consumption to GDP, investment to GDP and government expenditures to GDP ratios:

c_y = zeros(t);
x_y = zeros(t);
g_y = zeros(t);

c_y .= c_sim./y_sim;
x_y .= x_sim./y_sim;
g_y .= g_sim./y_sim;

Plot the ratios over time:

plot(1:t,c_y,label="Consumption to GDP",xlabel="Periods",title="GDP Decomposition in ratios")
plot!(1:t,x_y,label="Investment to GDP")
plot!(1:t,g_y,label="Gov. Exp. to GDP")

Lastly, recall from the firms FOCs that prices depend directly on TFP shocks. Let's plot a graph of interest rate, wages, and TFP shocks to see how those variables behave throughout the business cycle:

plot(1:t,z_sim,label="TFP shocks",xlabel="Periods",title="TFP, Interest Rates and Wages")
plot!(1:t,r_sim,label="Interest rates")
plot!(1:t,w_sim,label="Wages")