Solving and Simulating RBC Models with Mathematica

Diallo Ibrahima Amadou*
Clermont University, University of Auvergne, Centre d’Études et de Recherches sur le Développement International, CERDI-CNRS, 65, bd François Mitterrand, 63000 Clermont-Ferrand, France. 
Contact: zavren@gmail.com
September 2012

Abstract

This notebook illustrates how to solve and simulate a Real Business Cycle (RBC) model using Mathematica. It starts from finding the first order conditions, then derives the steady-state, log-linearizes and solves the RBC model. After solving the model, the notebook first computes the impulse response functions and then simulates the model. It ends by calculating some second moments statistics of the RBC model. When we want to solve a RBC model, soon or later, we will need to use a software. This notebook demonstrates how to do this using Mathematica without having the need to use paper and pencil to derive the equilibrium conditions and log-linearize the model first, and then use a numerical software to perform the simulations. In this notebook, all these steps are done with Mathematica and fully explained in detail.

Keywords: real business cycle; steady-state; log-linearization; simulation; impulse response functions; second moments statistics; dynamic stochastic general equilibrium models

JEL Classification: C6, E2, O3

*All comments are welcome. All Errors and inaccuracies are mine.

1. Introduction

Real Business Cycle (RBC) models have become the basic workhorse of Dynamic Stochastic General Equilibrium Models (DSGE). DSGE models themselves are now being used both by researchers and professionals, including central bankers, to study and forecast short term economic fluctuations. The hoped contribution of the notebook is the illustration of how to solve and simulate RBC models using Mathematica. Before starting, it is important to inform the reader that this notebook is written with Mathematica 8.0 in a Windows 7 machine. Thus the notebook may or may not work well in future versions of Mathematica.

2. Defining the Model

We load the PlotLegends package. This package allows us to make legends for graphs further below.

<<PlotLegends`

We present a modified version of the Hansen (1985) RBC model. The social planner’s maximization problem can be written as:

Max (Subscript[E, 0])[(∑_(t=0)^∞)(β^t)((((Subscript[c, t]^(1-θ))-1)/(1-θ))-ψ ((Subscript[l, t]^(1+γ))/(1+γ)))]
Subject to:
(Subscript[c, t])+(Subscript[iv, t])=(Subscript[y, t])
       (Subscript[y, t])=(Subscript[z, t]) (Subscript[k, t]^α) (Subscript[l, t]^(1-α))
       (Subscript[k, t+1])=(1-δ) (Subscript[k, t])+(Subscript[iv, t])
       Log[(Subscript[z, t+1])]=(1-ρ) Log[zs]+ρ Log[(Subscript[z, t])]+(Subscript[ϵ, t+1])
       (Subscript[k, 0]) and (Subscript[z, 0]) are given

Where:
(Subscript[E, 0])[·] is the expectations operator
  β ∈ (0,1) is the discount factor
  (Subscript[c, t]) is the consumption
  (1/θ) measures the constant intertemporal elasticity of substitution in consumption. We have θ >0
  γ represents the inverse of the Frisch intertemporal elasticity of substitution in labour supply. We have γ≥0
  (Subscript[l, t]) is labor 
  ψ indicates the disutility weight given to labor in the utility function. We have ψ>0
  (Subscript[iv, t]) is investment. We denote investment by (Subscript[iv, t]) since I is a reserved word in Mathematica
  (Subscript[y, t]) is output
  (Subscript[k, t]) is capital at period t and (Subscript[k, t+1]) is capital at t+1
  δ ∈ (0,1) is the capital depreciation rate 
  (Subscript[z, t]) is the exogenous technology shock. zs is the steady-state level of (Subscript[z, t])
  ρ ∈ [0,1] is the autocorrelation parameter of technology
  (Subscript[ϵ, t+1]) the innovations to technology
  (Subscript[ϵ, t])(Overscript[∼, iid])N(0,(σ^2))
  t indicates the time.

The model presented here results to the Hansen (1985) model when θ and γ tend towards 1 and 0 respectively. Let’s illustrate this:
First we take the limit of the utility function when θ goes towards 1

ssrbcmod1=Limit[(((Subscript[c, t]^(1-θ))-1)/(1-θ))-ψ ((Subscript[l, t]^(1+γ))/(1+γ)),θ->1]

Then we take the limit of the result when γ approaches 0.

ssrbcmod2=Limit[ssrbcmod1,γ->0]

This last expression is the utility function used in the Hansen (1985) RBC model. The only difference is that (Subscript[l, t]) represents there leisure instead of labor supply. Also the minus sign becomes a positive sign since we are dealing there with leisure.

Next we illustrate how the budget constraint is obtained. We first enter the equation of motion of capital stock.

ssrbcmod3=(Subscript[k, t+1])==(1-δ) (Subscript[k, t])+(Subscript[iv, t])

We enter the identity stating that consumption plus investment is equal to output.

ssrbcmod4=(Subscript[c, t])+(Subscript[iv, t])==(Subscript[y, t])

We solve for investment in the previous identity.

ssrbcmod5=Solve[ssrbcmod4,(Subscript[iv, t])][[1,1]]

We substitute investment in the equation of motion of capital.

ssrbcmod6=ssrbcmod3/.ssrbcmod5

We replace output by its value. Hence the budget constraint.

ssrbcmod7=ssrbcmod6/. (Subscript[y, t])->(Subscript[z, t]) (Subscript[k, t]^α) (Subscript[l, t]^(1-α))

2.1. The Bellman Equation

The Bellman equation is
v[(Subscript[k, t]),(Subscript[z, t])]=(Underscript[Max, ((Subscript[c, t]),(Subscript[k, t+1]),(Subscript[l, t]))]){(((Subscript[c, t]^(1-θ))-1)/(1-θ))-ψ ((Subscript[l, t]^(1+γ))/(1+γ))+β (Subscript[EX, t])[v[(Subscript[k, t+1]),(Subscript[z, t+1])]]}
Subject to:
The economy resource constraint: (Subscript[k, 1+t])=-(Subscript[c, t])+(1-δ) (Subscript[k, t])+(Subscript[l, t]^(1-α)) (Subscript[k, t]^α) (Subscript[z, t])
       And technology:  Log[(Subscript[z, t+1])]=(1-ρ) Log[zs]+ρ Log[(Subscript[z, t])]+(Subscript[ϵ, t+1])
We denote the expectations operator as (Subscript[EX, t])[·]since E is a reserved word in Mathematica.

We enter the left hand side of the Bellman equation

ssrbcmod8=v[(Subscript[k, t]),(Subscript[z, t])]

We enter the right hand side of the Bellman equation without the Max operator.

ssrbcmod9=(((Subscript[c, t]^(1-θ))-1)/(1-θ))-ψ ((Subscript[l, t]^(1+γ))/(1+γ))+β (Subscript[EX, t])[v[(Subscript[k, t+1]),(Subscript[z, t+1])]]

We enter the technology equation

ssrbcmod10=Log[(Subscript[z, t+1])]==(1-ρ) Log[zs]+ρ Log[(Subscript[z, t])]+(Subscript[ϵ, t+1])

We enter the right hand side of the Bellman equation without the expectations operator.

ssrbcmod11=(((Subscript[c, t]^(1-θ))-1)/(1-θ))-ψ ((Subscript[l, t]^(1+γ))/(1+γ))+β v[(Subscript[k, t+1]),(Subscript[z, t+1])]

We isolate consumption in the resource constraint.

ssrbcmod12=ssrbcmod7[[1]]-(Subscript[k, 1+t])==ssrbcmod7[[2]]-(Subscript[k, 1+t])

ssrbcmod13=ssrbcmod12[[1]]+(Subscript[c, t])==ssrbcmod12[[2]]+(Subscript[c, t])

We substitute consumption by its value in the right hand side of the Bellman Equation.

ssrbcmod14=ssrbcmod11/.ssrbcmod13[[1]]->ssrbcmod13[[2]]

We take the derivative of this last expression with respect to (Subscript[k, t+1]) and set it equal to 0.

ssrbcmod15=D[ssrbcmod14,(Subscript[k, t+1])]==0

We take the derivative of the right hand side of the Bellman Equation with respect to (Subscript[l, t]) and set it equal to 0.

ssrbcmod16=D[ssrbcmod14,(Subscript[l, t])]==0

The Benveniste-Scheinkman envelope theorem condition is given by

ssrbcmod17=D[ssrbcmod8,(Subscript[k, t])]==D[ssrbcmod14,(Subscript[k, t])]

We forward t by one period in the preceding equation.

ssrbcmod18=ssrbcmod17/.t->t+1

We take expections into account in equation ssrbcmod15

ssrbcmod19=ssrbcmod15/.Derivative[1,0][v][(Subscript[k, 1+t]),(Subscript[z, 1+t])]->(Subscript[EX, t])[Derivative[1,0][v][(Subscript[k, 1+t]),(Subscript[z, 1+t])]]

We substitute Derivative[1,0][v][(Subscript[k, 1+t]),(Subscript[z, 1+t])] by its value in the previous equation.

ssrbcmod20=ssrbcmod19/.ssrbcmod18[[1]]->ssrbcmod18[[2]]

Remember that consumption is given by

ssrbcmod21=ssrbcmod13

We forward t by one period in the preceding equation.

ssrbcmod22=ssrbcmod21/.t->t+1

We replace the right hand sides of ssrbcmod21 and ssrbcmod22 by (Subscript[c, t]) and (Subscript[c, t+1]) respectively in equation ssrbcmod20.

ssrbcmod23=ssrbcmod20/.{ssrbcmod21[[2]]->ssrbcmod21[[1]],ssrbcmod22[[2]]->ssrbcmod22[[1]]}

This last equation is known as the Euler equation. It is also called Lucas asset pricing equation. We rewrite it without the expectations operator. But we keep in mind that there is always an expectation operator for this equation.

ssrbcmod24=ssrbcmod23/.(Subscript[EX, t])[(Subscript[c, 1+t]^(-θ)) (1-δ+α (Subscript[k, 1+t]^(-1+α)) (Subscript[l, 1+t]^(1-α)) (Subscript[z, 1+t]))]->(Subscript[c, 1+t]^(-θ)) (1-δ+α (Subscript[k, 1+t]^(-1+α)) (Subscript[l, 1+t]^(1-α)) (Subscript[z, 1+t]))

Now, we replace (Subscript[c, t]) by its value in equation ssrbcmod16

ssrbcmod25=ssrbcmod16/.ssrbcmod21[[2]]->ssrbcmod21[[1]]

We solve for (Subscript[c, t]) in the previous equation

ssrbcmod26=Solve[ssrbcmod25,(Subscript[c, t])]

We transform this expression to an equation for convenience

ssrbcmod27=ssrbcmod26[[1,1,1]]==ssrbcmod26[[1,1,2]]

This last equation represents the labor-leisure choice equation. We expand the powers in this equation to get:

ssrbcmod28=PowerExpand[ssrbcmod27]

2.2. The Steady-State Equilibrium

Now we derive the steady-state equilibrium for the main variables. We start with the Euler equation.

ssrbcmod29=ssrbcmod23

We define steady-state variables by adding an s to the existing variables. For example the steady-state consumption is noted as cs for consumption star.

ssrbcmod30=ssrbcmod29/.{(Subscript[c, t])->cs,(Subscript[c, t+1])->cs,(Subscript[k, t+1])->ks,(Subscript[l, t+1])->ls,(Subscript[z, t+1])->zs}

By the definition of steady-steady, all variables are constant. Hence their expectations are themselves. Thus we have.

ssrbcmod31=ssrbcmod30/.(Subscript[EX, t])[(cs^(-θ)) (1+(ks^(-1+α)) (ls^(1-α)) zs α-δ)]->(cs^(-θ)) (1+(ks^(-1+α)) (ls^(1-α)) zs α-δ)

Also at the steady-state, there is no uncertainty, thus zs = 1.

ssrbcmod32=ssrbcmod31/.zs->1

We divide this last equation by (cs^(-θ))

ssrbcmod33=Expand[((ssrbcmod32[[1]])/(cs^(-θ)))==((ssrbcmod32[[2]])/(cs^(-θ)))]

We define the following ratio

ssrbcmod34=Solve[(ks/ls)==kls,ks][[1,1]]

We substitute ks by its value in ssrbcmod33

ssrbcmod35=ssrbcmod33/.ssrbcmod34//PowerExpand

We solve for kls in the previous equation.

ssrbcmod36=Solve[ssrbcmod35,kls][[1,1]]

We transform this result to an equation.

ssrbcmod37=ssrbcmod36[[1]]==ssrbcmod36[[2]]

Now we turn to the labor-leisure choice equation.

ssrbcmod38=ssrbcmod28

We define the steady-state variables for this equation.

ssrbcmod39=ssrbcmod38/.{(Subscript[c, t])->cs,(Subscript[l, t])->ls,(Subscript[k, t])->ks,(Subscript[z, t])->zs}

We set zs = 1

ssrbcmod40=ssrbcmod39/.zs->1

We substitute ks by its value

ssrbcmod41=ssrbcmod40/.ssrbcmod34//PowerExpand

We replace kls by its value.

ssrbcmod42=ssrbcmod41/.ssrbcmod36

We turn to the resource constraint.

ssrbcmod43=ssrbcmod7

ssrbcmod44=ssrbcmod43/.{(Subscript[k, t+1])->ks,(Subscript[c, t])->cs,(Subscript[k, t])->ks,(Subscript[l, t])->ls,(Subscript[z, t])->1}

We divide by ks.

ssrbcmod45=Map[(1/ks) #&,ssrbcmod44]//Expand

We substitute ks by its value.

ssrbcmod46=ssrbcmod45/.ssrbcmod34//PowerExpand

We replace kls and cs by their values.

ssrbcmod47=ssrbcmod46/.{ssrbcmod36,ssrbcmod42[[1]]->ssrbcmod42[[2]]}

We solve for ls in the previous equation.

ssrbcmod48=Solve[ssrbcmod47,ls]//FullSimplify

We change this to an equation.

ssrbcmod49=ssrbcmod48[[1,1,1]]==ssrbcmod48[[1,1,2]]

We find the value of cs.

ssrbcmod50=ssrbcmod42/.ssrbcmod49[[1]]->ssrbcmod49[[2]]

ssrbcmod51=FullSimplify[ssrbcmod50]

We recall the value of ks.

ssrbcmod52=ssrbcmod34

We write this as an equation.

ssrbcmod53=ssrbcmod52[[1]]==ssrbcmod52[[2]]

We then find the value of ks.

ssrbcmod54=ssrbcmod53/.{ssrbcmod36,ssrbcmod49[[1]]->ssrbcmod49[[2]]}

ssrbcmod55=FullSimplify[ssrbcmod54]

Next we find the value of ys.

ssrbcmod56=(Subscript[y, t])==(Subscript[z, t]) (Subscript[k, t]^α) (Subscript[l, t]^(1-α))

ssrbcmod57=ssrbcmod56/.{(Subscript[y, t])->ys,(Subscript[z, t])->zs,(Subscript[k, t])->ks,(Subscript[l, t])->ls}

We replace ks, ls and zs by their respective values.

ssrbcmod58=ssrbcmod57/.{ssrbcmod55[[1]]->ssrbcmod55[[2]],ssrbcmod49[[1]]->ssrbcmod49[[2]],zs->1}

ssrbcmod59=FullSimplify[ssrbcmod58]

We recall the expression of (Subscript[iv, t]).

ssrbcmod60=ssrbcmod5

We transform this expression to an equation.

ssrbcmod61=ssrbcmod60[[1]]==ssrbcmod60[[2]]

We find the value of ivs.

ssrbcmod62=ssrbcmod61/.{(Subscript[iv, t])->ivs,(Subscript[c, t])->cs,(Subscript[y, t])->ys}

We replace cs and ys by their respective values.

ssrbcmod63=ssrbcmod62/.{ssrbcmod51[[1]]->ssrbcmod51[[2]],ssrbcmod59[[1]]->ssrbcmod59[[2]]}

3. Log-Linearization

Log-linearization is a mean of transforming the dynamic nonlinear equations to dynamic linear equations. This makes the calculations and interpretations easier. To do the log-linearization, we will write a variable (Subscript[x, t]) as xs (e^(Subscript[xh, t])), where xs is the steady-state value of (Subscript[x, t]) and (Subscript[xh, t])=Log[((Subscript[x, t])/xs)]. This last expression is approximately the percentage deviation of (Subscript[x, t]) around the steady-state xs. In fact:

ssrbcmod64=(Subscript[xh, t])==Series[Log[((Subscript[x, t])/xs)],{(Subscript[x, t]),xs,1}]//Normal

The h in (Subscript[xh, t]) is put for x hat. At the steady-state, we can verify that (Subscript[xh, t])=((-xs+(Subscript[x, t]))/xs)=0. In fact at the steady-state (Subscript[xh, t])=((-xs+xs)/xs)=0. We can check that the expression xs (e^(Subscript[xh, t])) gives the original variable (Subscript[x, t]) by calculating:

ssrbcmod65=xs(E^(Subscript[xh, t]))/.(Subscript[xh, t])->Log[((Subscript[x, t])/xs)]

More log-linearization techniques can be found in Uhlig (1997) and Zietz (2008). Now let's start the log-linearization of the key equilibrium equations.

3.1. The Labor-Leisure Choice Equation

Let's write a function that performs multivariate linear Taylor series approximation. The function takes three arguments. The first represents the mathematical function we want to linearize. The second is the list of variables with respect to which we want to linearize and the third argument is the list of points around which we want to take the linearization.

mulLinearTaylor[f_,x_,a_]:=Module[{fa,gradf,gradfa},fa=f//.Thread[x->a];gradf=D[f,{x}];gradfa=gradf//.Thread[x->a];fa+gradfa.(x-a)]

We recall the labor-leisure choice equation.

ssrbcmod66=ssrbcmod28

We first replace each variable in this equation by its equivalent transformed expression.

ssrbcmod67=ssrbcmod66/.{(Subscript[c, t])->cs (E^(Subscript[ch, t])),(Subscript[l, t])->ls (E^(Subscript[lh, t])),(Subscript[k, t])->ks (E^(Subscript[kh, t])),(Subscript[z, t])->zs (E^(Subscript[zh, t]))}

We take the first order Taylor series approximation of the left hand side of this equation with respect to (Subscript[ch, t]) around 0.

ssrbcmod68=Series[ssrbcmod67[[1]],{(Subscript[ch, t]),0,1}]//Normal

We take the first order multivariate Taylor series approximation of the right hand side of equation ssrbcmod67 with respect to (Subscript[lh, t]), (Subscript[zh, t]) and (Subscript[kh, t]) respectively around 0, 0 and 0.

ssrbcmod69=mulLinearTaylor[ssrbcmod67[[2]],{(Subscript[lh, t]),(Subscript[zh, t]),(Subscript[kh, t])},{0,0,0}]

We equalize the two expressions ssrbcmod68 and ssrbcmod69 to form an equation.

ssrbcmod70=ssrbcmod68==ssrbcmod69

Now, we form the steady-state equation of the labor-leisure choice equation.

ssrbcmod71=ssrbcmod66/.{(Subscript[c, t])->cs,(Subscript[l, t])->ls,(Subscript[k, t])->ks,(Subscript[z, t])->zs}

We substract equation ssrbcmod71 from equation ssrbcmod70.

ssrbcmod72=ssrbcmod70[[1]]-ssrbcmod71[[1]]==ssrbcmod70[[2]]-ssrbcmod71[[2]]

We divide equation ssrbcmod72 by equation ssrbcmod71.

ssrbcmod73=((ssrbcmod72[[1]])/(ssrbcmod71[[1]]))==Apply[Plus,Table[FullSimplify[((ssrbcmod72[[2,i]])/(ssrbcmod71[[2]]))],{i,1,Length[ssrbcmod72[[2]]]}]]

This last equality represents the log-linearized equation of the labor-leisure choice equation. In traditional economic notation, this equation can be rewritten as:
(Subscript[(OverHat[c]), t])=((α (Subscript[(OverHat[k]), t]))/θ)-(((α+γ) (Subscript[(OverHat[l]), t]))/θ)+((Subscript[(OverHat[z]), t])/θ)
Where a hat over a variable indicates percentage deviation of the variable around the steady-state.

3.2. The Euler Equation or Lucas Asset Pricing Equation

We recall the Euler equation without the expectations operator. We will put back this operator when we finish the log-linearization process.

ssrbcmod74=ssrbcmod24

We isolate (Subscript[c, t]).

ssrbcmod75=ssrbcmod74[[1]]-β (Subscript[c, 1+t]^(-θ)) (1-δ+α (Subscript[k, 1+t]^(-1+α)) (Subscript[l, 1+t]^(1-α)) (Subscript[z, 1+t]))==ssrbcmod74[[2]]-β (Subscript[c, 1+t]^(-θ)) (1-δ+α (Subscript[k, 1+t]^(-1+α)) (Subscript[l, 1+t]^(1-α)) (Subscript[z, 1+t]))

We multiply this last equality by -1.

ssrbcmod76=Map[-1 #&,ssrbcmod75]

We replace each variable in this equation by its equivalent transformed expression.

ssrbcmod77=ssrbcmod76/.{(Subscript[c, t])->cs (E^(Subscript[ch, t])),(Subscript[c, t+1])->cs (E^(Subscript[ch, t+1])),(Subscript[k, t+1])->ks (E^(Subscript[kh, t+1])),(Subscript[l, t+1])->ls (E^(Subscript[lh, t+1])),(Subscript[z, t+1])->zs (E^(Subscript[zh, t+1]))}

We take the first order Taylor series approximation of the left hand side of this equation with respect to (Subscript[ch, t]) around 0.

ssrbcmod78=Series[ssrbcmod77[[1]],{(Subscript[ch, t]),0,1}]//Normal

We take the first order multivariate Taylor series approximation of the right hand side of equation ssrbcmod77 with respect to (Subscript[ch, 1+t]), (Subscript[zh, 1+t]), (Subscript[kh, 1+t]) and (Subscript[lh, 1+t]) respectively around 0, 0, 0 and 0.

ssrbcmod79=mulLinearTaylor[ssrbcmod77[[2]],{(Subscript[ch, 1+t]),(Subscript[zh, 1+t]),(Subscript[kh, 1+t]),(Subscript[lh, 1+t])},{0,0,0,0}]

We equalize the two expressions ssrbcmod78 and ssrbcmod79 to form an equation.

ssrbcmod80=ssrbcmod78==ssrbcmod79

Now, we form the steady-state equation of the Euler equation.

ssrbcmod81=ssrbcmod76/.{(Subscript[c, t])->cs,(Subscript[c, t+1])->cs,(Subscript[k, t+1])->ks,(Subscript[l, t+1])->ls,(Subscript[z, t+1])->zs}

We substract equation ssrbcmod81 from equation ssrbcmod80 side by side.

ssrbcmod82=ssrbcmod80[[1]]-ssrbcmod81[[1]]==ssrbcmod80[[2]]-ssrbcmod81[[2]]

We divide both sides of the previous equation by -(cs^(-θ)).

ssrbcmod83=((ssrbcmod82[[1]])/(-(cs^(-θ))))==Apply[Plus,Table[FullSimplify[((ssrbcmod82[[2,i]])/(-(cs^(-θ))))],{i,1,Length[ssrbcmod82[[2]]]}]]

We divide both sides of the steady-state equation by (cs^(-θ)).

ssrbcmod84=Map[(1/(cs^(-θ))) #&,ssrbcmod81]

We rewrite this equation as:

ssrbcmod85=ssrbcmod84[[2]]==ssrbcmod84[[1]]

This last equation implies that the log-linearized Euler equation is :

ssrbcmod86=ssrbcmod83/.ssrbcmod85[[1]]->ssrbcmod85[[2]]

We set zs = 1 in this last equation.

ssrbcmod87=ssrbcmod86/.zs->1

This last equation represents the log-linearized Euler equation without the expectations operator. This equality implies that the log-linearized Euler equation with the expectations operator is:

ssrbcmod88=ssrbcmod87/.{(Subscript[ch, 1+t])->(Subscript[EX, t])[(Subscript[ch, 1+t])],(Subscript[kh, 1+t])->(Subscript[EX, t])[(Subscript[kh, 1+t])],(Subscript[lh, 1+t])->(Subscript[EX, t])[(Subscript[lh, 1+t])],(Subscript[zh, 1+t])->(Subscript[EX, t])[(Subscript[zh, 1+t])]}

In traditional economic notation, this equation can be rewritten as:
θ (Subscript[(OverHat[c]), t])=θ (Subscript[E, t])((Subscript[(OverHat[c]), t+1]))-(ks^(-1+α)) (ls^(1-α)) (-1+α) α β (Subscript[E, t])((Subscript[(OverHat[k]), t+1]))+(ks^(-1+α)) (ls^(1-α)) (-1+α) α β (Subscript[E, t])((Subscript[(OverHat[l]), t+1]))-(ks^(-1+α)) (ls^(1-α)) α β (Subscript[E, t])((Subscript[(OverHat[z]), t+1]))
Where a hat over a variable indicates percentage deviation of the variable around the steady-state; (Subscript[E, t])(·) is the expectations operator and an s near a variable designates the steady-state value of that variable.

3.3. The Resource Constraint

We recall the resource constraint.

ssrbcmod89=ssrbcmod7

We replace each variable in this equation by its equivalent transformed expression.

ssrbcmod90=ssrbcmod89/.{(Subscript[k, t+1])->ks (E^(Subscript[kh, t+1])),(Subscript[c, t])->cs (E^(Subscript[ch, t])),(Subscript[k, t])->ks (E^(Subscript[kh, t])),(Subscript[l, t])->ls (E^(Subscript[lh, t])),(Subscript[z, t])->zs (E^(Subscript[zh, t]))}

We take the first order Taylor series approximation of the left hand side of this equation with respect to (Subscript[kh, 1+t]) around 0.

ssrbcmod91=Series[ssrbcmod90[[1]],{(Subscript[kh, 1+t]),0,1}]//Normal

We take the first order multivariate Taylor series approximation of the right hand side of equation ssrbcmod90 with respect to (Subscript[ch, t]), (Subscript[zh, t]), (Subscript[kh, t]) and (Subscript[lh, t]) respectively around 0, 0, 0 and 0.

ssrbcmod92=mulLinearTaylor[ssrbcmod90[[2]],{(Subscript[ch, t]),(Subscript[zh, t]),(Subscript[kh, t]),(Subscript[lh, t])},{0,0,0,0}]

We equalize the two expressions ssrbcmod91 and ssrbcmod92 to form an equation.

ssrbcmod93=ssrbcmod91==ssrbcmod92

Now, we form the steady-state equation of the Resource constraint.

ssrbcmod94=ssrbcmod89/.{(Subscript[k, t+1])->ks,(Subscript[c, t])->cs,(Subscript[k, t])->ks,(Subscript[l, t])->ls,(Subscript[z, t])->zs}

We substract equation ssrbcmod94 from equation ssrbcmod93 side by side.

ssrbcmod95=ssrbcmod93[[1]]-ssrbcmod94[[1]]==ssrbcmod93[[2]]-ssrbcmod94[[2]]

We divide both sides of this equation by ks.

ssrbcmod96=((ssrbcmod95[[1]])/ks)==Apply[Plus,Table[FullSimplify[((ssrbcmod95[[2,i]])/ks)],{i,1,Length[ssrbcmod95[[2]]]}]]

We set zs = 1 in this last equation.

ssrbcmod97=ssrbcmod96/.zs->1

This last equation represents the log-linearized resource constraint equation. In traditional economic notation, this equation can be rewritten as:
(Subscript[(OverHat[k]), t+1])=-((cs (Subscript[(OverHat[c]), t]))/ks)+(1+(ks^(-1+α)) (ls^(1-α)) α-δ) (Subscript[(OverHat[k]), t])-(ks^(-1+α)) (ls^(1-α)) (-1+α) (Subscript[(OverHat[l]), t])+(ks^(-1+α)) (ls^(1-α)) (Subscript[(OverHat[z]), t])
Where a hat over a variable indicates percentage deviation of the variable around the steady-state and an s near a variable designates the steady-state value of that variable.

3.4. The Technology Equation

We recall the technology equation.

ssrbcmod98=ssrbcmod10

We replace each variable in this equation by its equivalent transformed expression.

ssrbcmod99=ssrbcmod98/.{(Subscript[z, 1+t])->zs (E^(Subscript[zh, t+1])),(Subscript[z, t])->zs (E^(Subscript[zh, t]))}

We take the first order Taylor series approximation of the left hand side of this equation with respect to (Subscript[zh, 1+t]) around 0.

ssrbcmod100=Series[ssrbcmod99[[1]],{(Subscript[zh, 1+t]),0,1}]//Normal

We take the first order Taylor series approximation of the right hand side of equation ssrbcmod99 with respect to (Subscript[zh, t]) around 0.

ssrbcmod101=Series[ssrbcmod99[[2]],{(Subscript[zh, t]),0,1}]//Normal

We equalize the two expressions ssrbcmod100 and ssrbcmod101 to form an equation.

ssrbcmod102=ssrbcmod100==ssrbcmod101

Now, we form the steady-state equation of the Technology equation.

ssrbcmod103=ssrbcmod98/.{(Subscript[z, 1+t])->zs,(Subscript[z, t])->zs,(Subscript[ϵ, 1+t])->0}

We substract equation ssrbcmod103 from equation ssrbcmod102 side by side.

ssrbcmod104=ssrbcmod102[[1]]-ssrbcmod103[[1]]==ssrbcmod102[[2]]-ssrbcmod103[[2]]//Simplify

We rewrite this result as:

ssrbcmod105=ssrbcmod104[[2]]==ssrbcmod104[[1]]

This last equation represents the log-linearized technology equation. In traditional economic notation, this equation can be rewritten as:

(Subscript[(OverHat[z]), t+1])=ρ (Subscript[(OverHat[z]), t])+(Subscript[ϵ, 1+t])
Where a hat over a variable indicates percentage deviation of the variable around the steady-state. Taking expectations of both sides of the technology equation, we get:

ssrbcmod106=(Subscript[EX, t])[(Subscript[zh, 1+t])]==ρ (Subscript[zh, t])

We rewrite this last equation without the expectations operator but we keep in mind that there is always an expectations operator with this equation written in this way.

ssrbcmod107=(Subscript[zh, 1+t])==ρ (Subscript[zh, t])

3.5. The Output Equation

We recall the equation of output.

ssrbcmod108=ssrbcmod56

We replace each variable in this equation by its equivalent transformed expression.

ssrbcmod109=ssrbcmod108/.{(Subscript[y, t])->ys (E^(Subscript[yh, t])),(Subscript[k, t])->ks (E^(Subscript[kh, t])),(Subscript[l, t])->ls (E^(Subscript[lh, t])),(Subscript[z, t])->zs (E^(Subscript[zh, t]))}

We take the first order Taylor series approximation of the left hand side of this equation with respect to (Subscript[yh, t]) around 0.

ssrbcmod110=Series[ssrbcmod109[[1]],{(Subscript[yh, t]),0,1}]//Normal

We take the first order multivariate Taylor series approximation of the right hand side of equation ssrbcmod109 with respect to (Subscript[zh, t]), (Subscript[kh, t]) and (Subscript[lh, t]) respectively around 0, 0 and 0.

ssrbcmod111=mulLinearTaylor[ssrbcmod109[[2]],{(Subscript[zh, t]),(Subscript[kh, t]),(Subscript[lh, t])},{0,0,0}]

We equalize the two expressions ssrbcmod110 and ssrbcmod111 to form an equation.

ssrbcmod112=ssrbcmod110==ssrbcmod111

Now, we form the steady-state equation of the Output equation.

ssrbcmod113=ssrbcmod108/.{(Subscript[y, t])->ys,(Subscript[k, t])->ks,(Subscript[l, t])->ls,(Subscript[z, t])->zs}

We substract equation ssrbcmod113 from equation ssrbcmod112 side by side.

ssrbcmod114=ssrbcmod112[[1]]-ssrbcmod113[[1]]==ssrbcmod112[[2]]-ssrbcmod113[[2]]

We divide equation ssrbcmod114 by equation ssrbcmod113 side by side.

ssrbcmod115=((ssrbcmod114[[1]])/(ssrbcmod113[[1]]))==((ssrbcmod114[[2]])/(ssrbcmod113[[2]]))//Expand

We collect (Subscript[lh, t]).

ssrbcmod116=Collect[ssrbcmod115,(Subscript[lh, t])]

This last equation represents the log-linearized equation of output. In traditional economic notation, this equation can be rewritten as:

(Subscript[(OverHat[y]), t])=α (Subscript[(OverHat[k]), t])+(1-α) (Subscript[(OverHat[l]), t])+(Subscript[(OverHat[z]), t])
Where a hat over a variable indicates percentage deviation of the variable around the steady-state.

3.6. The Identity Consumption plus Investment is Equal to Output

We recall the expression of Investment.

ssrbcmod117=ssrbcmod61

We replace each variable in this equation by its equivalent transformed expression.

ssrbcmod118=ssrbcmod117/.{(Subscript[iv, t])->ivs (E^(Subscript[ivh, t])),(Subscript[c, t])->cs (E^(Subscript[ch, t])),(Subscript[y, t])->ys (E^(Subscript[yh, t]))}

We take the first order Taylor series approximation of the left hand side of this equation with respect to (Subscript[ivh, t]) around 0.

ssrbcmod119=Series[ssrbcmod118[[1]],{(Subscript[ivh, t]),0,1}]//Normal

We take the first order multivariate Taylor series approximation of the right hand side of equation ssrbcmod118 with respect to (Subscript[ch, t]) and (Subscript[yh, t]) respectively around 0 and 0.

ssrbcmod120=mulLinearTaylor[ssrbcmod118[[2]],{(Subscript[ch, t]),(Subscript[yh, t])},{0,0}]

We equalize the two expressions ssrbcmod119 and ssrbcmod120 to form an equation.

ssrbcmod121=ssrbcmod119==ssrbcmod120

Now, we form the steady-state equation of the Identity.

ssrbcmod122=ssrbcmod117/.{(Subscript[iv, t])->ivs,(Subscript[c, t])->cs,(Subscript[y, t])->ys}

We substract equation ssrbcmod122 from equation ssrbcmod121 side by side.

ssrbcmod123=ssrbcmod121[[1]]-ssrbcmod122[[1]]==ssrbcmod121[[2]]-ssrbcmod122[[2]]

We divide both sides of this equation by ivs.

ssrbcmod124=Map[(1/ivs) #&,ssrbcmod123]//Expand

This last equation represents the log-linearized equation of the Identity Consumption plus Investment is Equal to Output. In traditional economic notation, this equation can be rewritten as:

(Subscript[(OverHat[i]), t])=-((cs (Subscript[(OverHat[c]), t]))/ivs)+((ys (Subscript[(OverHat[y]), t]))/ivs)
Where a hat over a variable indicates percentage deviation of the variable around the steady-state and an s near a variable designates the steady-state value of that variable.

4. Solving the Model

We start by defining the matrix form representation of the model.

4.1. Matrix Form Representation

Following Hartley, Hoover and Salyer (1998), we write the model as a rational expectations equilibrium system of difference equation using the following equations: the labor-leisure choice equation, the Euler equation, the resource constraint and the technology equation. First, we define the vectors (Subscript[u, t])=((Subscript[(OverHat[c]), t]),(Subscript[(OverHat[k]), t]),(Subscript[(OverHat[l]), t]),(Subscript[(OverHat[z]), t])) and (Subscript[u, t+1])=((Subscript[(OverHat[c]), t+1]),(Subscript[(OverHat[k]), t+1]),(Subscript[(OverHat[l]), t+1]),(Subscript[(OverHat[z]), t+1])). With this notation, the system of difference equation can be written as A (Subscript[u, t])=B (Subscript[E, t])((Subscript[u, t+1])), where A and B are matrices to be found and (Subscript[E, t])(·) is the rational expectations operator. 
We enter the vectors (Subscript[u, t]) and (Subscript[u, t+1]).

ssrbcmod125=(Subscript[u, t])=={(Subscript[ch, t]),(Subscript[kh, t]),(Subscript[lh, t]),(Subscript[zh, t])}

ssrbcmod126=ssrbcmod125/.t->t+1

We remind the four log-linearized main equations of our model.

ssrbcmod127=ssrbcmod73
ssrbcmod128=ssrbcmod87
ssrbcmod129=ssrbcmod97
ssrbcmod130=ssrbcmod107

We cast these equations in the form of A (Subscript[u, t])=B (Subscript[E, t])((Subscript[u, t+1])) by putting all t variables on the left hand side and all t + 1 variables on the right hand side.

ssrbcmod131=ssrbcmod127[[1]]-ssrbcmod127[[2]]==0
ssrbcmod132=ssrbcmod128
ssrbcmod133=ssrbcmod129[[2]]==ssrbcmod129[[1]]
ssrbcmod134=ssrbcmod130[[2]]==ssrbcmod130[[1]]

Now, we start extracting the A matrix.

ssrbcmod135=CoefficientArrays[{ssrbcmod131,ssrbcmod132,ssrbcmod133,ssrbcmod134},ssrbcmod125[[2]]]//Normal

We take the A matrix.

ssrbcmod136=A==ssrbcmod135[[2]]

We make a copy of the previous expression.

ssrbcmod137=ssrbcmod136

We display this result in matrix form.

ssrbcmod138=MatrixForm[ssrbcmod137]//TraditionalForm

Now, we start extracting the B matrix.

ssrbcmod139=CoefficientArrays[{ssrbcmod131,ssrbcmod132,ssrbcmod133,ssrbcmod134},ssrbcmod126[[2]]]//Normal

We extract the second term of this result.

ssrbcmod140=ssrbcmod139[[2]]

We multiply this last result by -1 to obtain the B matrix.

ssrbcmod141=B==Map[-1 #&,ssrbcmod140,{2}]

We make a copy of the previous expression.

ssrbcmod142=ssrbcmod141

We display this result in matrix form.

ssrbcmod143=MatrixForm[ssrbcmod142]//TraditionalForm

4.2. Defining the Parameters (Calibration)

In this section, we define the values of the parameters (calibration). We take these values from the literature. θ = 3 is from Barro and Sala-i-Martin (2004); γ = 0.33 is from Mandelman and Zanetti (2008); β = 0.99,  α = 0.36, δ = 0.025, ψ = 3, (σ^2)=0.007 (so σ = 0.083666) and  ρ = 0.95  are from Hartley, Hoover and Salyer (1998). With these values of the parameters we can find the values of the steady-state variables, all equations and, the matrices A and B. 
We start by finding the numerical steady-state value of labor. We remind the equation of the steady-state value of labor.

ssrbcmod144=ssrbcmod49

Then we have:

ssrbcmod145=ssrbcmod144/.{α->0.36,β->0.99,δ->0.025,θ->3,ψ->3,γ->0.33}

We remind the steady-state value of consumption.

ssrbcmod146=ssrbcmod51

We find the numerical steady-state value of consumption.

ssrbcmod147=ssrbcmod146/.{α->0.36,β->0.99,δ->0.025,θ->3,ψ->3,γ->0.33}

We remind the steady-state value of capital.

ssrbcmod148=ssrbcmod55

We find the numerical steady-state value of capital.

ssrbcmod149=ssrbcmod148/.{α->0.36,β->0.99,δ->0.025,θ->3,ψ->3,γ->0.33}

We remind the steady-state value of output.

ssrbcmod150=ssrbcmod59

We find the numerical steady-state value of output.

ssrbcmod151=ssrbcmod150/.{α->0.36,β->0.99,δ->0.025,θ->3,ψ->3,γ->0.33}

We remind the steady-state value of investment.

ssrbcmod152=ssrbcmod63

We find the numerical steady-state value of investment.

ssrbcmod153=ssrbcmod152/.{α->0.36,β->0.99,δ->0.025,θ->3,ψ->3,γ->0.33}

An alternative way to find the numerical steady-state values is to enter the values of the parameters with an equal sign at the begining of this section and recall all the equations we want. But we avoid this method since we might still need symbolic values of the equations in the remaining of the notebook. Now we find the numerical values of the matrices A and B. We first remind the value of matrix A.

ssrbcmod154=ssrbcmod136

Then we find the numerical value of matrix A.

ssrbcmod155=ssrbcmod154/.{α->0.36,δ->0.025,θ->3,γ->0.33,ρ->0.95,cs->ssrbcmod147[[2]],ks->ssrbcmod149[[2]],ls->ssrbcmod145[[2]]}//N

We make a copy of matrix A.

ssrbcmod156=ssrbcmod155

We display this last result in matrix form.

ssrbcmod157=MatrixForm[ssrbcmod156]//TraditionalForm

We remind the value of matrix B.

ssrbcmod158=ssrbcmod141

Then we find the numerical value of matrix B.

ssrbcmod159=ssrbcmod158/.{α->0.36,θ->3,β->0.99,ks->ssrbcmod149[[2]],ls->ssrbcmod145[[2]]}//N

We make a copy of matrix B.

ssrbcmod160=ssrbcmod159

We display this last result in matrix form.

ssrbcmod161=MatrixForm[ssrbcmod160]//TraditionalForm

To solve the model, we will follow the solution method proposed by Hartley, Hoover and Salyer (1998). Remember that we wrote the model as A (Subscript[u, t])=B (Subscript[E, t])((Subscript[u, t+1])). Premultiply both sides of this system by (A^(-1)). We get (Subscript[u, t])=(A^(-1)) B (Subscript[E, t])((Subscript[u, t+1])). Then we decompose the matrix (A^(-1)) B as: (A^(-1)) B = Q Λ (Q^(-1)). Where Q is a matrix which has as columns the eigenvectors of the matrix (A^(-1)) B and Λ is a diagonal matrix which has as diagonal elements the eigenvalues of the matrix (A^(-1)) B. Replacing the matrix (A^(-1)) B by its value, we get: (Subscript[u, t])=Q Λ(Q^(-1)) (Subscript[E, t])((Subscript[u, t+1])). We premultiply both sides of this equation by (Q^(-1)). We then have:
 (Q^(-1)) (Subscript[u, t])=Λ(Q^(-1)) (Subscript[E, t])((Subscript[u, t+1])). We rewrite this as:  (Q^(-1)) (Subscript[u, t])=Λ (Subscript[E, t])((Q^(-1))(Subscript[u, t+1])). We set (Subscript[d, t])=(Q^(-1))(Subscript[u, t]) and (Subscript[d, t+1]) = (Q^(-1))(Subscript[u, t+1]). Hence we have: 
 (Subscript[d, t]) = Λ (Subscript[E, t])((Subscript[d, t+1])). Since Λ is a diagonal matrix, we can rewrite this last equation as independent equations having the following form: (Subscript[d, i,t])= (Subscript[λ, i]) (Subscript[E, t])((Subscript[d, i,t+1])). Where i indicates the equations and (Subscript[λ, i]) is an eigenvalue of the matrix (A^(-1)) B. It is important to note that i = 1..4 here since (Q^(-1)) is 4 by 4 matrix in this model. Performing recursive forward substitution in the equations (Subscript[d, i,t])= (Subscript[λ, i]) (Subscript[E, t])((Subscript[d, i,t+1])) for T periods yields: (Subscript[d, i,t])=(Subscript[λ, i]^T) (Subscript[E, t])((Subscript[d, i,t+T])). This last equation implies that the stable rational expectations solution to our model is related with the eigenvalue that has an absolute value less than unity. This result comes from the fact that taking the limit of the expression (Subscript[λ, i]^T) (Subscript[E, t])((Subscript[d, i,t+T])) as T goes toward infinity, gives 0 only if |(Subscript[λ, i])| is less than one. This implies that (Subscript[d, i,t]) = 0 for that particular eigenvalue. This last expression is the solution we are looking for.
Now let’s find the (A^(-1)) B matrix.

ssrbcmod162=Inverse[ssrbcmod155[[2]]].ssrbcmod159[[2]]

Then we find the eigenvalues and eigenvectors of matrix (A^(-1)) B.

ssrbcmod163=Eigensystem[ssrbcmod162]

We extract the eigenvalues from the previous result.

ssrbcmod164=ssrbcmod163[[1]]

We form the Λ matrix.

ssrbcmod165=Λ==DiagonalMatrix[ssrbcmod164]

We make a copy of the previous result.

ssrbcmod166=ssrbcmod165

We display this last result in matrix form.

ssrbcmod167=MatrixForm[ssrbcmod166]//TraditionalForm

This last result shows us that the third eigenvalue equal to 0.946224 has an absolute value less than one. This result implies that the stable rational expectations solution to our model is related to this eigenvalue.

Now we extract the eigenvectors from expression ssrbcmod163.

ssrbcmod168=ssrbcmod163[[2]]

We form the matrix (Q^(-1)).

ssrbcmod169=Inverse[ssrbcmod168]//Transpose

We make a copy of the previous result.

ssrbcmod170=ssrbcmod169

We display (Q^(-1)) in matrix form.

ssrbcmod171=MatrixForm[ssrbcmod170]//TraditionalForm

As we explained previously, a stable rational expectations solution is related to the eigenvalue less than one in absolute value which is here the third eigenvalue. This implies that the third row of the matrix (Q^(-1))gives the linear restriction we are looking for. Hence we extract the third row of matrix (Q^(-1)).

ssrbcmod172=ssrbcmod169[[3]]

Remember that the rational expectations solution is given by the expression (Subscript[d, i,t]) = 0. Also recall that (Subscript[d, t])=(Q^(-1))(Subscript[u, t]). Hence we form the equation (Subscript[d, i,t]) = 0, which is our rational expectations solution.

ssrbcmod173=ssrbcmod172.ssrbcmod125[[2]]==0

We isolate (Subscript[ch, t]) in this last expression.

ssrbcmod174=Solve[ssrbcmod173,(Subscript[ch, t])]//Expand

We transform this expression to an equation. This is our rational expectations solution.

ssrbcmod175=ssrbcmod174[[1,1,1]]==ssrbcmod174[[1,1,2]]

In addition to this last equation, the other equations describing the equilibrium conditions of our model are: the economy resource constraint, the labor-leisure choice equation, the equation of output, the identity consumption plus investment is equal to output and the technology equation. In what follows, we find the numerical values of these equations.

We recall the economy resource constraint.

ssrbcmod176=ssrbcmod97

We calculate the numerical value of the resource constraint.

ssrbcmod177=ssrbcmod176/.{α->0.36,δ->0.025,cs->ssrbcmod147[[2]],ks->ssrbcmod149[[2]],ls->ssrbcmod145[[2]]}

We recall the labor-leisure choice equation.

ssrbcmod178=ssrbcmod73

We isolate (Subscript[lh, t]) in this last equation.

ssrbcmod179=Solve[ssrbcmod178,(Subscript[lh, t])]//Expand

We transform this expression to an equation.

ssrbcmod180=ssrbcmod179[[1,1,1]]==ssrbcmod179[[1,1,2]]

We find the numerical value of the labor-leisure choice equation.

ssrbcmod181=ssrbcmod180/.{θ->3,α->0.36,γ->0.33}

We recall the equation of output.

ssrbcmod182=ssrbcmod116

We find the numerical value of the equation of output.

ssrbcmod183=ssrbcmod182/.α->0.36

We recall the identity consumption plus investment is equal to output.

ssrbcmod184=ssrbcmod124

We find the numerical value of the identity consumption plus investment is equal to output.

ssrbcmod185=ssrbcmod184/.{cs->ssrbcmod147[[2]],ivs->ssrbcmod153[[2]],ys->ssrbcmod151[[2]]}

We recall the technology equation.

ssrbcmod186=ssrbcmod105

We lag t by one period in this last equation.

ssrbcmod187=ssrbcmod186/.t->t-1

We find the numerical value of the technology equation.

ssrbcmod188=ssrbcmod187/.ρ->0.95

We express the variables (Subscript[ch, t]),(Subscript[lh, t]),(Subscript[yh, t]),(Subscript[kh, 1+t]) and (Subscript[ivh, t]) in function of (Subscript[kh, t]) and (Subscript[zh, t]) only.

ssrbcmod189=Solve[{ssrbcmod175,ssrbcmod177,ssrbcmod181,ssrbcmod183,ssrbcmod185},{(Subscript[ch, t]),(Subscript[lh, t]),(Subscript[yh, t]),(Subscript[kh, 1+t]),(Subscript[ivh, t])}]

We rearrange this last expression as:

ssrbcmod190=Chop[ssrbcmod189]//Flatten

Now we transform each expression of this preceding result to an equation. We begin with the equation of consumption.

ssrbcmod191=ssrbcmod190[[1,1]]==ssrbcmod190[[1,2]]

The equation of labor.

ssrbcmod192=ssrbcmod190[[2,1]]==ssrbcmod190[[2,2]]

The equation of output.

ssrbcmod193=ssrbcmod190[[3,1]]==ssrbcmod190[[3,2]]

The equation of next period capital ((Subscript[kh, 1+t])).

ssrbcmod194=ssrbcmod190[[4,1]]==ssrbcmod190[[4,2]]

The equation of investment.

ssrbcmod195=ssrbcmod190[[5,1]]==ssrbcmod190[[5,2]]

Now we write functions that numerically evaluate each of these preceding equations. We start with the technology equation. We first remind the technology equation.

ssrbcmod196=ssrbcmod188

The function that numerically evaluates the technology equation is given in the next line. This function takes three arguments: the first is the autocorrelation parameter, the second represents a list containing the numerical values of the innovations to technology and the third is the total time length we want for the technology series.

fzh[ρ_,e_,T_]:=Module[{res1,zh,t,T1},T1=T+1;zh=Table[0,{i,T1}];res1=Table[0,{i,T1}];
For[t=2,t≤T1,t++,
zh[[t]]=ρ zh[[t-1]]+e[[t]];
res1[[t]]=zh[[t]];
];
res1]

We recall the equation of next period capital stock ((Subscript[kh, 1+t])).

ssrbcmod197=ssrbcmod194

Now we write a function that calculates the series for the capital stock ((Subscript[kh, t])). This function takes four arguments: the first two arguments represent the numerical coefficients of the equation of next period capital, the third argument is the total time length we want for the capital stock series and the fourth is the technology series.

fkh[a_,b_,T_,zh_]:=Module[{res1,res2,t,T1,kh},T1=T+1;kh=Table[0,{i,T1}];res1=Table[0,{i,T1}];
For[t=1,t≤(T1-1),t++,
kh[[t+1]]=a kh[[t]]+b zh[[t]];
res1[[t+1]]=kh[[t+1]];
];
res2=res1;Most[res2]]

We remind the equation of consumption.

ssrbcmod198=ssrbcmod191

We write a function that calculates the consumption series. This function takes four arguments: the first two are respectively the coefficients of capital and technology in the equation of consumption, the third is the capital series and the fourth the technology series.

fch[ac_,bc_,kh_,zh_]:=ac kh+bc zh

We remind the equation of labor.

ssrbcmod199=ssrbcmod192

We write a function that calculates the labor series. This function takes four arguments: the first two are respectively the coefficients of capital and technology in the equation of labor, the third is the capital series and the fourth the technology series.

flh[al_,bl_,kh_,zh_]:=al kh+bl zh

We remind the equation of output.

ssrbcmod200=ssrbcmod193

We write a function that calculates the output series. This function takes four arguments: the first two are respectively the coefficients of capital and technology in the equation of output, the third is the capital series and the fourth the technology series.

fyh[ay_,by_,kh_,zh_]:=ay kh+by zh

We remind the equation of investment.

ssrbcmod201=ssrbcmod195

We write a function that calculates the investment series. This function takes four arguments: the first two are respectively the coefficients of capital and technology in the equation of investment, the third is the capital series and the fourth the technology series.

fivh[aiv_,biv_,kh_,zh_]:=aiv kh+biv zh

5. Studying the Economy

In this section we study the economy formed by our theoretical model.

5.1. Impulse Response Functions

To generate the impulse response functions, we first create the innovations to technology. We shock these innovations by one standard deviation and plot the result.

ϵp1=Table[0,{i,100}];

ϵp2=Insert[ϵp1,0.083666,2];

ssrbcmod202=ListLinePlot[ϵp2,FrameLabel->{"Years after shock","(Subscript[ϵ, t])"},PlotLabel->"Shock to the innovations of Technology",Frame->True,GridLines->Automatic,GridLinesStyle->Dashed]

Now we remind the equation of technology and calculate the impulse response function of technology.

ssrbcmod196

ssrbcmod203=fzh[0.95,ϵp2,99];

We plot the impulse response function of technology. The graph shows that technology starts from its steady-state value of 0 and increases in response to the shock of one standard deviation, and then gradually decreases towards its steady-state value.

ssrbcmod204=ListLinePlot[ssrbcmod203,FrameLabel->{"Years after shock","% deviation from steady-state"},PlotLabel->"Impulse Response of Technology, (Subscript[(OverHat[z]), t])",Frame->True,GridLines->Automatic,GridLinesStyle->Dashed,PlotStyle->{Thickness[0.007],Hue[0.65]}]

We remind the equation of next period capital stock ((Subscript[kh, t+1])) and calculate the impulse response function of capital stock ((Subscript[kh, t])).

ssrbcmod197

ssrbcmod205=fkh[ssrbcmod197[[2,1,1]],ssrbcmod197[[2,2,1]],100,ssrbcmod203];

We plot the impulse response function of capital stock ((Subscript[kh, t])). We see that the graph of the impulse response function of capital stock is hump-shaped.

ssrbcmod206=ListLinePlot[ssrbcmod205,FrameLabel->{"Years after shock","% deviation from steady-state"},PlotLabel->"Impulse Response of Capital Stock, (Subscript[(OverHat[k]), t])",Frame->True,GridLines->Automatic,GridLinesStyle->Dashed,PlotStyle->{Thickness[0.007],Hue[0.1]}]

We remind the equation of consumption.

ssrbcmod198

We calculate the impulse response function of consumption.

ssrbcmod207=fch[ssrbcmod198[[2,1,1]],ssrbcmod198[[2,2,1]],ssrbcmod205,ssrbcmod203];

We plot the impulse response function of consumption. Like the capital stock, we see that the graph of the impulse response function of consumption is hump-shaped.

ssrbcmod208=ListLinePlot[ssrbcmod207,FrameLabel->{"Years after shock","% deviation from steady-state"},PlotLabel->"Impulse Response of Consumption, (Subscript[(OverHat[c]), t])",Frame->True,GridLines->Automatic,GridLinesStyle->Dashed,PlotStyle->{Thickness[0.007],Hue[0.3]}]

We remind the equation of labor.

ssrbcmod192

We calculate the impulse response function of labor.

ssrbcmod209=flh[ssrbcmod192[[2,1,1]],ssrbcmod192[[2,2,1]],ssrbcmod205,ssrbcmod203];

We plot the impulse response function of labor. We see that labor first augments and then decreases gradually, and increases again progressively to its steady-state level.

ssrbcmod210=ListLinePlot[ssrbcmod209,FrameLabel->{"Years after shock","% deviation from steady-state"},PlotLabel->"Impulse Response of Labor, (Subscript[(OverHat[l]), t])",Frame->True,GridLines->Automatic,GridLinesStyle->Dashed,PlotStyle->{Thickness[0.007],Hue[0.5]}]

We remind the equation of output.

ssrbcmod200

We calculate the impulse response function of output.

ssrbcmod211=fyh[ssrbcmod200[[2,1,1]],ssrbcmod200[[2,2,1]],ssrbcmod205,ssrbcmod203];

We plot the impulse response function of output. As technology, we observe that output increases in response to the shock and then gradually falls to its steady-state level.

ssrbcmod212=ListLinePlot[ssrbcmod211,FrameLabel->{"Years after shock","% deviation from steady-state"},PlotLabel->"Impulse Response of Output, (Subscript[(OverHat[y]), t])",Frame->True,GridLines->Automatic,GridLinesStyle->Dashed,PlotStyle->{Thickness[0.007],Hue[0.8]}]

We remind the equation of investment.

ssrbcmod201

We calculate the impulse response function of investment.

ssrbcmod213=fivh[ssrbcmod201[[2,1,1]],ssrbcmod201[[2,2,1]],ssrbcmod205,ssrbcmod203];

We plot the impulse response function of investment. As output, we see that investment rises in response to the shock and then progressively falls towards its steady-state level.

ssrbcmod214=ListLinePlot[ssrbcmod213,FrameLabel->{"Years after shock","% deviation from steady-state"},PlotLabel->"Impulse Response of Investment, (Subscript[(OverHat[i]), t])",Frame->True,GridLines->Automatic,GridLinesStyle->Dashed,PlotStyle->{Thickness[0.007],Hue[0.99]},PlotRange->{{0,100},{-0.1,0.35}}]

Now we combine all these preceding figures into a single graph.

ssrbcmod215=ListLinePlot[{ssrbcmod203,ssrbcmod205,ssrbcmod207,ssrbcmod209,ssrbcmod211,ssrbcmod213},FrameLabel->{"Years after shock","% deviation from steady-state"},PlotLabel->"Impulse Response Functions of the Variables",Frame->True,GridLines->Automatic,GridLinesStyle->Dashed,PlotLegend->{"Technology","Capital Stock","Consumption","Labor","Output","Investment"}, LegendPosition->{1.1,-0.4},ImageSize->550,PlotStyle->{{Thickness[0.007],Hue[0.65]},{Thickness[0.007],Hue[0.1]},{Thickness[0.007],Hue[0.3]},{Thickness[0.007],Hue[0.5]},{Thickness[0.007],Hue[0.8]},{Thickness[0.007],Hue[0.99]}},PlotRange->{{0,100},{-0.05,0.35}}]

Now we present the graphs side by side.

ssrbcmod216=GraphicsGrid[{{ssrbcmod204,ssrbcmod206},{ssrbcmod208,ssrbcmod210},{ssrbcmod212,ssrbcmod214}},ImageSize->600]

5.2. Stochastic Simulation

We set the random seed to 2000 in order to allow the reproduction of our simulations.

SeedRandom[2000]

We generate 5000 normal random numbers with mean 0 and standard deviation 0.083666. This series represents the innovations to technology.

ssrbcmod217=RandomVariate[NormalDistribution[0,0.083666],5000];

Now we plot only the last 150 values of the previous series.

ssrbcmod218=ListLinePlot[Take[ssrbcmod217,-150],FrameLabel->{"Time",(Subscript[ϵ, t])},PlotLabel->"Innovations to technology",Frame->True,GridLines->Automatic,GridLinesStyle->Dashed]

Next we present the Hodrick-Prescott (1997) filter used to filter the data we are going to study. To find the Hodrick-Prescott (HP) filter, we minimize the following expression:

(Underscript[Min, (Subscript[((Subscript[g, t])), t=1]^T)])(∑_(t=1)^T)(Subscript[c, t]^2)+λ (∑_(t=3)^T)(([((Subscript[g, t])-(Subscript[g, t-1]))-((Subscript[g, t-1])-(Subscript[g, t-2]))])^2)

Subject to:

  (Subscript[c, t])=(Subscript[y, t])-(Subscript[g, t])

Where (Subscript[y, t]) is the series we want to smooth, (Subscript[g, t]) is the trend component, (Subscript[c, t]) is the cyclical component and λ is the smoothing parameter. 
Danthine and Girardin (1989) show that the solution to this optimization problem is given by:

(Subscript[g, t])=(((I+λ K' K))^(-1)) (Subscript[y, t])

With I an identity matrix of dimension T and K is a special matrix of dimension (T-2)×T having the following form:

(({{({{({{({{1}, {0}})}, {0}})}, {⋮}})}, {0}})({{({{({{({{-2}, {1}})}, {0}})}, {⋮}})}, {0}})({{({{({{({{1}, {-2}})}, {1}})}, {⋮}})}, {0}})({{({{({{({{0}, {1}})}, {-2}})}, {⋮}})}, {0}})({{({{({{({{0}, {0}})}, {1}})}, {⋮}})}, {0}})({{({{({{({{⋯}, {⋯}})}, {⋯}})}, {⋯}})}, {⋯}})({{({{({{({{0}, {0}})}, {0}})}, {⋮}})}, {1}})({{({{({{({{0}, {0}})}, {0}})}, {⋮}})}, {-2}})({{({{({{({{0}, {0}})}, {0}})}, {⋮}})}, {1}}))

Having found (Subscript[g, t]), the the cyclical component is calculated as  (Subscript[c, t])=(Subscript[y, t])-(Subscript[g, t]). Since the parameters of our model are calibrated for quarterly data, we will use the value λ = 1600. Now let us write a function that computes the HP filter. The function takes two arguments: the first is the data we want to smooth and the second is the value of the parameter λ. The function then returns the cyclical component and the trend component respectively.

hPfilterz[data_List,λ_?Positive]:=Module[{totime,trend,cycle,KKMAT,IDM,MATA},totime=Length[data];KKMAT=SparseArray[
    {{r_, r_} -> 1.,{r_, s_} /; s == r + 1 -> -2.,{r_, s_} /; s == r + 2 -> 1.}, {(totime - 2), totime},0.];IDM=IdentityMatrix[totime];MATA=IDM+λ Transpose[KKMAT].KKMAT ;trend = LinearSolve[MATA, data,Method->"Automatic"];cycle=data-trend;{cycle,trend}]

We generate the simulated time series of technology.

ssrbcmod219=fzh[0.95,ssrbcmod217,4999];

We take the last 150 values of the previous series.

ssrbcmod220=Take[ssrbcmod219,-150];

We apply the HP-filter to the preceding result.

ssrbcmod221=hPfilterz[ssrbcmod220,1600];

We extract the cyclical component of the technology series.

ssrbcmod222=ssrbcmod221[[1]];

We plot the HP-filtered time path of technology.

ssrbcmod223=ListLinePlot[ssrbcmod222,FrameLabel->{"Time","% deviation from steady-state"},PlotLabel->"Time path of Technology (HP-filtered), (Subscript[(OverHat[z]), t])",Frame->True,GridLines->Automatic,GridLinesStyle->Dashed,PlotStyle->{Thickness[0.004],Hue[0.65]}]

We generate the simulated time series of capital stock.

ssrbcmod224=fkh[ssrbcmod197[[2,1,1]],ssrbcmod197[[2,2,1]],5000,ssrbcmod219];

We take the last 150 values of the previous series.

ssrbcmod225=Take[ssrbcmod224,-150];

We apply the HP-filter to the preceding result.

ssrbcmod226=hPfilterz[ssrbcmod225,1600];

We extract the cyclical component of the capital stock series.

ssrbcmod227=ssrbcmod226[[1]];

We plot the HP-filtered time path of capital stock.

ssrbcmod228=ListLinePlot[ssrbcmod227,FrameLabel->{"Time","% deviation from steady-state"},PlotLabel->"Time path of Capital Stock (HP-filtered), (Subscript[(OverHat[k]), t])",Frame->True,GridLines->Automatic,GridLinesStyle->Dashed,PlotStyle->{Thickness[0.004],Hue[0.1]}]

We generate the simulated time series of consumption.

ssrbcmod229=fch[ssrbcmod198[[2,1,1]],ssrbcmod198[[2,2,1]],ssrbcmod224,ssrbcmod219];

We take the last 150 values of the previous series.

ssrbcmod230=Take[ssrbcmod229,-150];

We apply the HP-filter to the preceding result.

ssrbcmod231=hPfilterz[ssrbcmod230,1600];

We extract the cyclical component of the consumption series.

ssrbcmod232=ssrbcmod231[[1]];

We plot the HP-filtered time path of consumption.

ssrbcmod233=ListLinePlot[ssrbcmod232,FrameLabel->{"Time","% deviation from steady-state"},PlotLabel->"Time path of Consumption (HP-filtered), (Subscript[(OverHat[c]), t])",Frame->True,GridLines->Automatic,GridLinesStyle->Dashed,PlotStyle->{Thickness[0.004],Hue[0.3]}]

We generate the simulated time series of labor.

ssrbcmod234=flh[ssrbcmod192[[2,1,1]],ssrbcmod192[[2,2,1]],ssrbcmod224,ssrbcmod219];

We take the last 150 values of the previous series.

ssrbcmod235=Take[ssrbcmod234,-150];

We apply the HP-filter to the preceding result.

ssrbcmod236=hPfilterz[ssrbcmod235,1600];

We extract the cyclical component of the labor series.

ssrbcmod237=ssrbcmod236[[1]];

We plot the HP-filtered time path of labor.

ssrbcmod238=ListLinePlot[ssrbcmod237,FrameLabel->{"Time","% deviation from steady-state"},PlotLabel->"Time path of Labor (HP-filtered), (Subscript[(OverHat[l]), t])",Frame->True,GridLines->Automatic,GridLinesStyle->Dashed,PlotStyle->{Thickness[0.004],Hue[0.5]}]

We generate the simulated time series of output.

ssrbcmod239=fyh[ssrbcmod200[[2,1,1]],ssrbcmod200[[2,2,1]],ssrbcmod224,ssrbcmod219];

We take the last 150 values of the previous series.

ssrbcmod240=Take[ssrbcmod239,-150];

We apply the HP-filter to the preceding result.

ssrbcmod241=hPfilterz[ssrbcmod240,1600];

We extract the cyclical component of the output series.

ssrbcmod242=ssrbcmod241[[1]];

We plot the HP-filtered time path of output.

ssrbcmod243=ListLinePlot[ssrbcmod242,FrameLabel->{"Time","% deviation from steady-state"},PlotLabel->"Time path of Output (HP-filtered), (Subscript[(OverHat[y]), t])",Frame->True,GridLines->Automatic,GridLinesStyle->Dashed,PlotStyle->{Thickness[0.004],Hue[0.8]}]

We generate the simulated time series of investment.

ssrbcmod244=fivh[ssrbcmod201[[2,1,1]],ssrbcmod201[[2,2,1]],ssrbcmod224,ssrbcmod219];

We take the last 150 values of the previous series.

ssrbcmod245=Take[ssrbcmod244,-150];

We apply the HP-filter to the preceding result.

ssrbcmod246=hPfilterz[ssrbcmod245,1600];

We extract the cyclical component of the investment series.

ssrbcmod247=ssrbcmod246[[1]];

We plot the HP-filtered time path of investment.

ssrbcmod248=ListLinePlot[ssrbcmod247,FrameLabel->{"Time","% deviation from steady-state"},PlotLabel->"Time path of Investment (HP-filtered), (Subscript[(OverHat[i]), t])",Frame->True,GridLines->Automatic,GridLinesStyle->Dashed,PlotStyle->{Thickness[0.004],Hue[0.99]}]

Now we combine all these preceding figures into a single graph.

ssrbcmod249=ListLinePlot[{ssrbcmod222,ssrbcmod227,ssrbcmod232,ssrbcmod237,ssrbcmod242,ssrbcmod247},FrameLabel->{"Time","% deviation from steady-state"},PlotLabel->"Time path of the Variables (HP-filtered)",Frame->True,GridLines->Automatic,GridLinesStyle->Dashed,PlotLegend->{"Technology","Capital Stock","Consumption","Labor","Output","Investment"}, LegendPosition->{1.1,-0.4},ImageSize->550,PlotStyle->{{Thickness[0.007],Hue[0.65]},{Thickness[0.007],Hue[0.1]},{Thickness[0.007],Hue[0.3]},{Thickness[0.007],Hue[0.5]},{Thickness[0.007],Hue[0.8]},{Thickness[0.007],Hue[0.99]}},PlotRange->{{0,150},{-1.5,1.5}}]

The graph illustrates that investment is more volatile than output while capital stock and consumption are smoother than output.

Now we present the graphs side by side.

ssrbcmod250=GraphicsGrid[{{ssrbcmod223,ssrbcmod228},{ssrbcmod233,ssrbcmod238},{ssrbcmod243,ssrbcmod248}},ImageSize->600]

5.3. Second Moments Statistics

We start by writing a function that calculates the relative standard deviation. The function takes two arguments. Both arguments are lists or time series in our case. The result of the function is the standard deviation of the first series relative to the second.

relativesd[mvar_List,output_List]:=Module[{sdmvar,sdout},sdmvar=StandardDeviation[mvar];sdout=StandardDeviation[output];sdmvar/sdout]

We compute the standard deviation of each series relative output. To do this we automate the process by creating a function that takes as first argument a list of list containing the series and as second argument the output series. The function then gives as result the list of the relative standard deviations.

listrsdfunct[listrsd_,outrsd_]:=Module[{res1,NN},NN=Length[listrsd];res1=Table[0,{i,NN}];Do[res1[[i]]=relativesd[listrsd[[i]],outrsd],{i,1,NN}];res1]

We enter the series for which we want to calculate the relative standard deviation with output.

ssrbcmod251={ssrbcmod222,ssrbcmod227,ssrbcmod232,ssrbcmod237,ssrbcmod247};

We calculate the relative standard deviation of technology, capital stock, consumption, labor and investment respectively with output.

ssrbcmod252=listrsdfunct[ssrbcmod251,ssrbcmod242]

We create a function that compute the correlation of each series with output. The function takes as first argument a list of list containing the series and as second argument the output series. The result of the function is a list containing the correlation of each series with output.

listcorrfunct[listcorr_,outcorr_]:=Module[{rescorr,NNC},NNC=Length[listcorr];rescorr=Table[0,{i,NNC}];Do[rescorr[[i]]=Correlation[listcorr[[i]],outcorr],{i,1,NNC}];rescorr]

We compute the correlation of technology, capital stock, consumption, labor and investment respectively with output.

ssrbcmod253=listcorrfunct[ssrbcmod251,ssrbcmod242]

We display the relative standard deviation next to the correlation.

ssrbcmod254=Transpose[{ssrbcmod252,ssrbcmod253}]

Now we display the Second Moments Statistics of the model in table form.

ssrbcmod255=TableForm[ssrbcmod254,TableHeadings->{{"Technology","Capital Stock","Consumption","Labor", "Investment"},{"relative volatility","correlation"}},TableAlignments->Center]

From this table it appears that consumption, labor and capital stock are less volatile than output. Technology is almost as volatile as output. Investment is more volatile than output. The table also illustrates that technology, investment and consumption are highly procyclic, while capital stock is less correlated with output.

6. Conclusion

This notebook has demonstrated how to solve and simulate a Real Business Cycle (RBC) model using the software Mathematica. We have also shown how to find the first order conditions, derive the steady-state and log-linearize the model. A good extension of this notebook could be:

First, to automate some parts instead of doing it step by step.We have chosen this last method for pedagogic reasons.

Second, to use alternatives solution methods.

7. References

Barro, R. J. and X. Sala-i-Martin: 2004, Economic Growth, Second Edition, MIT Press.

Danthine, J. P. and M. Girardin: 1989, Business cycles in switzerland, a comparative study, European Economic Review, 33, 31-50.

Hansen, G. D.: 1985, Indivisible labor and the business cycle, Journal of Monetary Economics, 16(3), 309-327.

Hartley, J. E., K. D. Hoover and K. D. Salyer: 1998, A user’s guide to solving real business cycle models, In Real business cycles, A Reader, pp: 43-54, Routledge, London and New York.

Hodrick, R. J. and E. C. Prescott: 1997, Post-war U.S. business cycles: an empirical investigation, Journal of Money, Credit and Banking, 29(1), 1-16.

Mandelman, F. S. and F. Zanetti: 2008, Estimating general equilibrium models: an application with labour market frictions, Technical Handbook-No.1, Centre for Central Banking Studies, Bank of England.

Uhlig, H.: 1997, A toolkit for analyzing nonlinear dynamic stochastic models easily, Working Paper, Center, University of Tilburg.

Zietz, J.: 2008, A clarifying note on converting to log-deviations from the steady state, Economics Bulletin, 3(50), 1-15.