In [1]:
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import statsmodels.api as sm
%matplotlib inline
Consider the following business cycle model: \begin{align} Y_t & = A_t K_t^{\alpha} \tag{1}\\ C_t & = (1-s)Y_t \tag{2}\\ I_t & = K_{t+1} - ( 1- \delta) \tag{3}\\ Y_t & = C_t + I_t \tag{4} \end{align} where: \begin{align} \log A_{t+1} & = \rho \log A_t + \epsilon_t, \tag{5} \end{align} reflects exogenous fluctuation in TFP. The endogenous variables in the model are $K_t$, $Y_t$, $C_t$, $I_t$, and $A_t$ and $\epsilon_t$ is an exogenous white noise shock process with standard deviation $\sigma$. $K_t$ and $A_t$ are called state variables because their values in period $t$ affect the equilibrium of the model in period $t+1$.
In the (non-stochastic) steady state: \begin{align} \epsilon_t & = 0 \end{align} and \begin{align} K_t & = 0 \end{align} for all $t$. So we drop the $t$ subscripts and write the steady state solution solution to the model as: \begin{align} A & = 1\\ K & = \left(\frac{sA}{\delta}\right)^{\frac{1}{1-\alpha}}\\ Y & = AK^{\alpha}\\ I & = Y - C \end{align}
In [2]:
# Define parameters
s = 0.1
delta = 0.025
alpha = 0.35
# Compute the steady state values of the endogenous variables
Kss = (s/delta)**(1/(1-alpha))
Yss = Kss**alpha
Css = (1-s)*Yss
Iss = Yss - Css
print('Steady states:\n')
print('capital: ',round(Kss,5))
print('output: ',round(Yss,5))
print('consumption:',round(Css,5))
print('investment: ',round(Iss,5))
Now, you will simulate how the model behaves in the presence of a set of random TFP shocks. The simulation will run for $T+1$ periods from $t = 0,\ldots, T$. Suppose that $T = 100$.
Initialize an array for $\epsilon_t$ called eps that is equal to a $T\times 1$ array of normally distributed random variables with mean 0 and standard deviation $\sigma = 0.006$. Set the seed for the Numpy random number generator to 192. Plot $\epsilon_t$.
Initialize an array for $\log A_t$ called logA that is equal to a $(T+1)\times 1$ array of zeros. Set $\rho = 0.75$ and use the simulated values for $\epsilon_t$ to compute $\log A_1, \log A_2, \ldots, \log A_T$. Plot $\log A_t$.
Create a new variable called A that stores simulated values of $A_t$ (Note: $A_t = e^{\log A_t}$). Plot $A_t$.
Initialize an array for $K_t$ called K that is a $(T+1)\times 1$ array of zeros. Set the first value in the array equal to steady state capital. Then compute the subsequent values for $K_t$ using the computed values for $A_t$. Plot $K_t$.
Create variables called Y, C, and I that store simulated values for $Y_t$, $C_t$, and $I_t$.
Construct a $2\times2$ grid of subplots of the simulated paths of capital, output, consumption, and investment.
Compute the log deviation of each variable from its steady state ($(X_t - X_{ss})/X_{ss}$) and store the results in variables called: k_dev, y_dev, c_dev, and i_dev.
Construct a $2\times2$ grid of subplots of the impulse responses of capital, output, consumption, and investment to the technology shock with each variable expressed as a deviation from steady state.
Save the simulated data in a DataFrame called data_simulated with columns output, consumption, investment, and tfp.
In [3]:
# Step 1: simulate eps
T = 100
sigma = 0.006
np.random.seed(192)
eps = np.random.normal(scale=sigma,size=T)
plt.plot(np.arange(1,T+1),eps,lw=3,alpha = 0.7)
plt.title('$\epsilon_t$')
plt.grid()
In [4]:
# Step 2: simulate and plot log(TFP) logA
rho = 0.75
logA = np.zeros(T+1)
for t in range(T):
logA[t+1] = rho*logA[t] + eps[t]
plt.plot(logA,lw=3,alpha = 0.7)
plt.title('$\log A_t$')
plt.grid()
In [5]:
# Step 3: compute and plot TFP A
A = np.exp(logA)
plt.plot(A,lw=3,alpha = 0.7)
plt.title('$A_t$')
plt.grid()
In [6]:
# Step 4: Compulte and plot capital K
K = np.zeros(T+1)
K[0] = Kss
for t in range(T):
K[t+1] = s*A[t]*K[t]**alpha + (1-delta)*K[t]
plt.plot(K,lw=3,alpha = 0.7)
plt.title('$K_t$')
plt.grid()
In [7]:
# Step 5: Compute Y, C, and I
Y = A*K**alpha
C = (1-s)*Y
I = s*Y
In [8]:
# Step 6: Create a 2x2 plot of y, c, i, and k
fig = plt.figure(figsize=(12,8))
ax = fig.add_subplot(2,2,1)
ax.plot(Y,lw=3,alpha = 0.7)
ax.set_title('$Y_t$')
ax.grid()
ax = fig.add_subplot(2,2,2)
ax.plot(C,lw=3,alpha = 0.7)
ax.set_title('$C_t$')
ax.grid()
ax = fig.add_subplot(2,2,3)
ax.plot(I,lw=3,alpha = 0.7)
ax.set_title('$I_t$')
ax.grid()
ax = fig.add_subplot(2,2,4)
ax.plot(K,lw=3,alpha = 0.7)
ax.set_title('$K_t$')
ax.grid()
plt.tight_layout()
In [9]:
# Step 7: Compute y_dev, c_dev, i_dev, and k_dev to be the log deviations from steady state of the
# respective variables
y_dev = np.log(Y) - np.log(Yss)
c_dev = np.log(C) - np.log(Css)
i_dev = np.log(I) - np.log(Iss)
k_dev = np.log(K) - np.log(Kss)
In [10]:
# Step 8: Create a 2x2 plot of y_dev, c_dev, i_dev, and k_dev
fig = plt.figure(figsize=(12,8))
ax = fig.add_subplot(2,2,1)
ax.plot(100*y_dev,lw=3,alpha = 0.7)
ax.set_title('$\hat{y}_t$')
ax.grid()
ax.set_ylabel('% dev from steady state')
ax = fig.add_subplot(2,2,2)
ax.plot(100*c_dev,lw=3,alpha = 0.7)
ax.set_title('$\hat{c}_t$')
ax.grid()
ax = fig.add_subplot(2,2,3)
ax.plot(100*i_dev,lw=3,alpha = 0.7)
ax.set_title('$\hat{\imath}_t$')
ax.grid()
ax.set_ylabel('% dev from steady state')
ax = fig.add_subplot(2,2,4)
ax.plot(100*k_dev,lw=3,alpha = 0.7)
ax.set_title('$\hat{k}_t$')
ax.grid()
plt.tight_layout()
In [11]:
# Step 9: Save the simulated data in a DataFrame called data_simulated
data_simulated = pd.DataFrame({'output':y_dev,'consumption':c_dev,'investment':i_dev,'tfp':logA})
data_simulated.head()
Out[11]:
In [12]:
# Create a DataFrame with actual cyclical components of output, consumption, investment, and TFP
data_actual = pd.read_csv('http://www.briancjenkins.com/teaching/winter2017/econ129/data/Econ129_Rbc_Data.csv',index_col=0)
data_actual = pd.DataFrame({'output':np.log(data_actual.gdp/data_actual.gdp_trend),
'consumption':np.log(data_actual.consumption/data_actual.consumption_trend),
'investment':np.log(data_actual.investment/data_actual.investment_trend),
'tfp':np.log(data_actual.tfp/data_actual.tfp_trend)})
data_actual.head()
Out[12]:
In [13]:
# Compute the standard deviations of the actual business cycle data
print(data_actual.std())
In [14]:
# Compute the standard deviations of the simulated business cycle data
print(data_simulated.std())
In [15]:
# Compute the coeficients of correlation for the actual business cycle data
print(data_actual.corr())
In [16]:
# Compute the coeficients of correlation for the actual business cycle data
print(data_simulated.corr())