In [1]:
%matplotlib inline
import numpy as np
import sympy as sy
import pandas as pd
import statsmodels.api as sm
import statsmodels.formula.api as smf
import matplotlib.pyplot as plt
from numpy.linalg import eig, inv
from pandas_datareader.data import DataReader
sy.init_printing(use_latex='mathjax')
np.set_printoptions(precision=4, suppress=True)
In [2]:
from matplotlib import rcParams
# Restore old behavior of rounding default axis ranges
rcParams['axes.autolimit_mode'] = 'round_numbers'
rcParams['axes.xmargin'] = 0
rcParams['axes.ymargin'] = 0
# Adjust tick placement
rcParams['xtick.direction'] = 'in'
rcParams['ytick.direction'] = 'in'
rcParams['xtick.top'] = True
rcParams['ytick.right'] = True
# Disable legend frame
rcParams['legend.frameon'] = False
In [3]:
rec = DataReader('USREC', 'fred', start='1947-01', end='2020-01')
rec = rec.dropna()
pop = DataReader(['B230RC0Q173SBEA','CNP16OV'], 'fred', start='1947-01', end='2020-01')
pop = pop.dropna()
pop = pop.resample('QS').mean()
pop = pop.dropna()
pop.tail()
Out[3]:
In [4]:
fred = DataReader(['PRS85006053','GDP','PCND','PCESV','PCDG',
'FPI','HOANBS','COMPNFB','GDPDEF','TB3MS'],
'fred', start='1947-01', end='2020-01')
fred = fred.dropna()
fred = fred.resample('QS').mean()
fred = fred.dropna()
fred.head()
Out[4]:
In [5]:
fred.tail()
Out[5]:
In [6]:
pop_hp_cycle, pop_hp_trend = sm.tsa.filters.hpfilter(np.log(pop), lamb=10000)
pop_smooth = np.exp(pop_hp_trend)
fig, ax = plt.subplots()
pop['CNP16OV'].to_period('D').plot(ax=ax, lw=2)
pop_smooth['CNP16OV'].to_period('D').plot(ax=ax, lw=2)
plt.legend(['raw', 'filtered'])
plt.show()
fig, ax = plt.subplots()
pop['CNP16OV'].pct_change().to_period('D').plot(ax=ax, lw=2)
pop_smooth['CNP16OV'].pct_change().to_period('D').plot(ax=ax, lw=2)
plt.legend(['raw', 'filtered'], loc='upper left')
plt.show()
In [7]:
dta = fred[['GDP']]
dta.is_copy = False
dta['Output'] = np.log(fred['GDP']*10**9/fred['GDPDEF']*100
/(pop_smooth['CNP16OV']*10**3))
dta['Consumption'] = np.log((fred['PCND']+fred['PCESV'])*10**9/fred['GDPDEF']*100
/(pop_smooth['CNP16OV']*10**3))
dta['Investment'] = np.log((fred['PCDG']+fred['FPI'])*10**9/fred['GDPDEF']*100
/(pop_smooth['CNP16OV']*10**3))
dta['Capital'] = 0*dta['Output']
dta['Hours'] = np.log(fred['HOANBS']*100*fred['GDP']/np.mean(fred['GDP']['2010-01':'2010-10'])
/fred['PRS85006053']/pop_smooth['CNP16OV'])
dta['Wages'] = np.log(fred['COMPNFB']/fred['GDPDEF']*100)
dta['Interest Rate'] = (1+fred['TB3MS']/100)**(1/4)/(1+fred['GDPDEF'].pct_change())
dta['TFP'] = 0*dta['Output']
dta['Productivity'] = np.log(fred['PRS85006053']/fred['GDPDEF']*100)-dta['Hours']
dta['Price Level'] = np.log(fred['GDPDEF'])
dta = dta.drop('GDP', 1)
dta = dta.dropna()
dta.head()
Out[7]:
In [8]:
dta.tail()
Out[8]:
In [9]:
# Estimate Capital series using PIM
temp = fred[['GDP']]
temp.is_copy = False
temp['Inv'] = (fred['PCDG']+fred['FPI'])/fred['GDPDEF']*100/4
temp['LnInv'] = np.log(temp['Inv'])
temp['t'] = np.arange(len(temp['Inv']))
trend = smf.ols(formula='LnInv ~ t', data=temp).fit()
intercept, slope = trend.params
delta = 0.025
K = np.zeros(len(temp['LnInv']))
K_init = np.exp(intercept-slope)/(slope+delta)
K[0] = (1-delta)*K_init+temp['Inv'][0]
for i in range(1,len(temp['Inv'])):
K[i] = (1-delta)*K[i-1]+temp['Inv'][i]
temp['Cap'] = K
temp['Cap'].to_period('D').plot(lw=2)
plt.show()
(temp['Cap']/(fred['GDP']/fred['GDPDEF']*100)).to_period('D').plot(lw=2)
plt.show()
In [10]:
dta['Capital'] = np.log(temp['Cap']*10**9/(pop_smooth['CNP16OV']*10**3))
α = 1/3
dta['TFP'] = dta['Output']-α*dta['Capital']-(1-α)*dta['Hours']
dta['TFP'].to_period('D').plot(lw=2)
plt.show()
In [11]:
RGDP_pc = (fred['GDP']*10**9/fred['GDPDEF']*100/(pop_smooth['B230RC0Q173SBEA']*10**3)).dropna()
fig, ax = plt.subplots()
RGDP_pc.to_period('D').plot(ax=ax, lw=2, style='k-')
ax.set_ylim(10000, 60000)
ylim = ax.get_ylim()
ax.fill_between(rec.index, ylim[0], ylim[1], rec['USREC'], facecolor='lightgrey', edgecolor='lightgrey')
plt.title('US real per capita GDP (2009 dollars)')
plt.show()
In [12]:
fig, ax = plt.subplots()
np.log(RGDP_pc).to_period('D').plot(ax=ax, lw=2, style='k-')
ticks = [10000, 20000, 30000, 40000, 50000, 60000]
ax.set_yticks(np.log(ticks))
ax.set_yticklabels(ticks)
ax.set_ylim(np.log(ticks[0]), np.log(ticks[-1]))
ylim = ax.get_ylim()
ax.set_ylim(ylim)
ax.fill_between(rec.index, ylim[0], ylim[1], rec['USREC'], facecolor='lightgrey', edgecolor='lightgrey')
plt.title('US real per capita GDP, log scale')
plt.show()
In [13]:
fig, ax = plt.subplots()
(100*RGDP_pc.pct_change(4)).to_period('D').plot(ax=ax, lw=2, style='k-')
avg = np.mean(100*RGDP_pc.to_period('D').pct_change(4))
ax.set_ylim(-6, 12)
ylim = ax.get_ylim()
ax.hlines(avg, dta.index[0], dta.index[-1], color='r', linewidth=2)
ax.hlines(0, dta.index[0], dta.index[-1], linewidth=0.5)
ax.fill_between(rec.index, ylim[0], ylim[1], rec['USREC'], facecolor='lightgrey', edgecolor='lightgrey')
plt.title('US real per capita GDP, year-over-year change (%)')
plt.show()
print('Average growth rate (%) =', avg)
In [14]:
gdp_lin_cycle, gdp_lin_trend = sm.tsa.filters.hpfilter(np.log(RGDP_pc), lamb=1e9)
gdp_hp_cycle, gdp_hp_trend = sm.tsa.filters.hpfilter(np.log(RGDP_pc))
gdp_cf_cycle, gdp_cf_trend = sm.tsa.filters.cffilter(np.log(RGDP_pc))
In [15]:
fig, ax = plt.subplots()
np.log(RGDP_pc).to_period('D').plot(ax=ax, lw=2, style='k-')
gdp_lin_trend.plot(ax=ax, lw=2, style='r-')
ticks = [10000, 20000, 30000, 40000, 50000, 60000]
ax.set_yticks(np.log(ticks))
ax.set_yticklabels(ticks)
ax.set_ylim(np.log(ticks[0]), np.log(ticks[-1]))
ylim = ax.get_ylim()
ax.set_ylim(ylim)
ax.fill_between(rec.index, ylim[0], ylim[1], rec['USREC'], facecolor='lightgrey', edgecolor='lightgrey')
plt.title('US real per capita GDP, exponential trend')
plt.show()
In [16]:
fig, ax = plt.subplots()
(100*gdp_lin_cycle).to_period('D').plot(ax=ax, lw=2, style='k-')
ax.hlines(0, dta.index[0], dta.index[-1], linewidth=0.5)
ylim = ax.get_ylim()
ax.set_ylim(ylim)
ax.fill_between(rec.index, ylim[0], ylim[1], rec['USREC'], facecolor='lightgrey', edgecolor='lightgrey')
plt.title('US real per capita GDP, exponential trend residual (%)')
plt.show()
In [17]:
fig, ax = plt.subplots()
np.log(RGDP_pc).to_period('D').plot(ax=ax, lw=2, style='k-')
gdp_hp_trend.plot(ax=ax, lw=2, style='r-')
ticks = [10000, 20000, 30000, 40000, 50000, 60000]
ax.set_yticks(np.log(ticks))
ax.set_yticklabels(ticks)
ax.set_ylim(np.log(ticks[0]), np.log(ticks[-1]))
ylim = ax.get_ylim()
ax.set_ylim(ylim)
ax.fill_between(rec.index, ylim[0], ylim[1], rec['USREC'], facecolor='lightgrey', edgecolor='lightgrey')
plt.title('US real per capita GDP, Hodrick-Prescott trend')
plt.show()
In [18]:
fig, ax = plt.subplots()
(100*gdp_hp_cycle.to_period('D')).plot(ax=ax, lw=2, style='k-')
ax.hlines(0, dta.index[0], dta.index[-1], linewidth=0.5)
ax.set_ylim(-6, 6)
ylim = ax.get_ylim()
ax.fill_between(rec.index, ylim[0], ylim[1], rec['USREC'], facecolor='lightgrey', edgecolor='lightgrey')
plt.title('US real per capita GDP, Hodrick-Prescott residual (%)')
plt.show()
In [19]:
fig, ax = plt.subplots()
(100*gdp_cf_cycle).to_period('D').plot(ax=ax, lw=2, style='k-')
ax.hlines(0, dta.index[0], dta.index[-1], linewidth=0.5)
ylim = ax.get_ylim()
ax.fill_between(rec.index, ylim[0], ylim[1], rec['USREC'], facecolor='lightgrey', edgecolor='lightgrey')
plt.title('US real per capita GDP, Christiano-Fitzgerald residual (%)')
plt.show()
In [20]:
np.corrcoef(gdp_hp_cycle, gdp_cf_cycle)
Out[20]:
In [21]:
hp_cycles, hp_trend = sm.tsa.filters.hpfilter((100*dta).dropna())
cf_cycles, cf_trend = sm.tsa.filters.cffilter((100*dta).dropna())
In [22]:
print('Standard Deviations')
print(hp_cycles.std())
print('')
print('Autocorrelations')
a = list(dta.columns.values)
for i in range(len(a)):
print(dta.columns.values[i], ' \t\t', hp_cycles[dta.columns.values[i]].autocorr())
print('')
print('Correlations')
print(hp_cycles.corr(method='pearson'))
In [23]:
cf_cycles['Investment / 4'] = cf_cycles['Investment'] / 4
fig, ((ax1, ax2, ax3), (ax4, ax5, ax6),
(ax7, ax8, ax9)) = plt.subplots(3, 3, figsize=(20,15), sharex=False, sharey=False)
cf_cycles[['Output','Consumption']].to_period('D').plot(ax=ax1, style=['k','r'])
cf_cycles[['Output','Investment / 4']].to_period('D').plot(ax=ax2, style=['k','r'])
cf_cycles[['Output','Capital']].to_period('D').plot(ax=ax3, style=['k','r'])
cf_cycles[['Output','Hours']].to_period('D').plot(ax=ax4, style=['k','r'])
cf_cycles[['Output','Wages']].to_period('D').plot(ax=ax5, style=['k','r'])
cf_cycles[['Output','Interest Rate']].to_period('D').plot(ax=ax6, style=['k','r'])
cf_cycles[['Output','TFP']].to_period('D').plot(ax=ax7, style=['k','r'])
cf_cycles[['Output','Productivity']].to_period('D').plot(ax=ax8, style=['k','r'])
cf_cycles[['Output','Price Level']].to_period('D').plot(ax=ax9, style=['k','r'])
plt.savefig('US_CF.pdf', transparent=True, bbox_inches='tight', pad_inches=0.05)
plt.show()
In [24]:
fig, ax = plt.subplots()
cf_cycles[['Output','Consumption','Investment / 4','Hours']].to_period('D').plot(ax=ax)
ylim = ax.get_ylim()
ax.fill_between(rec.index, ylim[0], ylim[1], rec['USREC'], facecolor='lightgrey', edgecolor='lightgrey')
plt.legend(ncol=2, frameon=True)
plt.show()
In [25]:
dta['t'] = np.arange(len(dta['TFP']))
trend_TFP = smf.ols(formula='TFP ~ t', data=dta).fit()
dta['TFP_resid'] = trend_TFP.resid
intercept_TFP, slope_TFP = trend_TFP.params
print(trend_TFP.summary())
In [26]:
dta['TFP'].to_period('D').plot(lw=2)
(intercept_TFP + slope_TFP*dta['t']).plot(lw=2)
plt.show()
In [27]:
dta['TFP_resid'].to_period('D').plot(lw=2)
plt.show()
In [28]:
resid_TFP = smf.ols(formula='TFP_resid ~ TFP_resid.shift() -1', data=dta).fit()
print(resid_TFP.summary())
print('')
print('TFP residual autocorrelation =', resid_TFP.params[0])
print('TFP residual standard deviation =', resid_TFP.resid.std())
In [29]:
class Dynare(object):
def __init__(self, var, varexo, param_values, model, initval):
self.var = var
self.varexo = varexo
self.param_values = param_values
self.model = model
self.initval = initval
self.TranslateInputs()
return None
def TranslateInputs(self):
"""
Convert input strings into sympy objects
and translate timing
"""
# Endogenous variables
self.var_symbols = sy.symbols(self.var)
self.n = len(self.var_symbols)
self.current = ()
self.future = ()
self.past = ()
self.steadys = ()
for i in range(self.n):
self.current += (sy.symbols(str(self.var_symbols[i])+str('_{t}'))),
self.future += (sy.symbols(str(self.var_symbols[i])+str('_{t+1}'))),
self.past += (sy.symbols(str(self.var_symbols[i])+str('_{t-1}'))),
self.steadys += (sy.symbols(str(self.var_symbols[i])+str('_{ss}'))),
# Exogenous variables
self.varexo_symbols = sy.symbols(self.varexo),
self.p = len(self.varexo_symbols)
self.shocks = ()
for i in range(self.p):
self.shocks += (sy.symbols(str(self.varexo_symbols[i])+str('_{t}'))),
# Model equations
self.symbol_dict = {'betta':sy.symbols('beta'), 'gama':sy.symbols('gamma')}
self.timing_dict = {}
for i in range(self.n):
self.timing_dict[sy.sympify(str(self.var_symbols[i])+'(+1)')] = self.future[i]
self.timing_dict[sy.sympify(str(self.var_symbols[i])+'(-1)')] = self.past[i]
self.timing_dict[sy.sympify(str(self.var_symbols[i]))] = self.current[i]
for i in range(self.p):
self.timing_dict[sy.sympify(str(self.varexo_symbols[i]))] = self.shocks[i]
self.model_symbols = sy.sympify(self.model)
self.model_symbols = self.model_symbols.subs(self.timing_dict)
self.model_symbols = self.model_symbols.subs(self.symbol_dict)
try:
temp = len(self.model_symbols)
except:
self.model_symbols = self.model_symbols,
self.system = sy.Matrix(self.model_symbols)
return None
def SteadySystem(self):
"""
Removes lead/lag structure from the model to prepare for steady state calculation
"""
self.steady_vars = {}
for i in range(self.n):
self.steady_vars[self.current[i]] = self.steadys[i]
self.steady_vars[self.future[i]] = self.steadys[i]
self.steady_vars[self.past[i]] = self.steadys[i]
for i in range(self.p):
self.steady_vars[self.shocks[i]] = 0
self.steady_system = self.system.subs(self.steady_vars)
return self.steady_system
def SteadyValues(self):
"""
Numerically solves for the steady state of the system given initval
"""
self.SteadySystem()
try:
ss = sy.nsolve(self.steady_system.subs(self.param_values), self.steadys, self.initval)
except:
raise RuntimeError('Adjust initial values')
ss = ss.T.tolist()
self.steady_values = {}
for i in range(self.n):
self.steady_values[self.steadys[i]] = np.float(ss[0][i])
return self.steady_values
def steady(self):
"""
Prints out steady state values for the user
"""
self.SteadyValues()
print('\n' + 'STEADY-STATE RESULTS' + '\n')
for i in range(self.n):
print(str(self.var_symbols[i]), '\t%.4f' % self.steady_values[self.steadys[i]])
return None
def resid(self):
self.SteadyValues()
temp = self.steady_system.subs(self.param_values).subs(self.steady_values)
print('\n' + 'Residuals of the static equations' + '\n')
for i in range(self.n):
print('Equation number', i, ': %.4f' % temp[i])
return None
def TimeIteration(self):
"""
Solves the first-order approximation of the model using time iteration (thanks to Pontus Rendahl)
"""
self.SteadyValues()
self.A_symb = self.system.jacobian(self.past)
self.B_symb = self.system.jacobian(self.current)
self.C_symb = self.system.jacobian(self.future)
self.D_symb = self.system.jacobian(self.shocks)
self.A = np.array(self.A_symb.subs(self.steady_vars).subs(self.steady_values).subs(self.param_values)).astype(float)
self.B = np.array(self.B_symb.subs(self.steady_vars).subs(self.steady_values).subs(self.param_values)).astype(float)
self.C = np.array(self.C_symb.subs(self.steady_vars).subs(self.steady_values).subs(self.param_values)).astype(float)
self.D = np.array(self.D_symb.subs(self.steady_vars).subs(self.steady_values).subs(self.param_values)).astype(float)
self.metric = 1
self.F = np.zeros((self.n, self.n))
self.S = np.zeros((self.n, self.n))
# Add maxit to while loop?
while self.metric > 1e-13:
self.F = inv(self.B + self.C @ self.F) @ (-self.A)
self.S = inv(self.B + self.A @ self.S) @ (-self.C)
self.metric = np.max(np.max(np.abs(self.A + self.B @ self.F + self.C @ self.F @ self.F)))
self.Q = -inv(self.B + self.C @ self.F) @ self.D
# Need formal BK check?
if sum(eig(self.F)[0] > 1) != 0:
raise RuntimeError('Blanchard Kahn conditions are not satisfied: no stable equilibrium')
if sum(eig(self.S)[0] > 1) != 0:
raise RuntimeError('Blanchard Kahn conditions are not satisfied: indeterminacy')
return None
def SimulatedMoments(self, hp_filter=None, shocks_stderr=0.01, periods=10000):
self.TimeIteration()
x = np.zeros((self.n, periods))
É› = np.zeros((self.p, periods))
for i in range(self.p):
É›[i, :] = shocks_stderr * np.random.randn(periods)
for t in range(1, periods):
x[:, t] = self.F @ x[:, t-1] + self.Q @ É›[:, t]
print('SIMULATED MOMENTS')
print('')
print('VARIABLE \t STD. DEV.')
if hp_filter == None:
for i in range(self.n):
print(str(self.var_symbols[i]), np.std(x[i, :]))
else:
self.SteadyValues()
try:
for i in range(self.n):
hp_cycle, hp_trend = sm.tsa.filters.hpfilter(100*np.log(x[i, :]+self.steady_values[self.steadys[i]]),
lamb=hp_filter)
print(str(self.var_symbols[i]), '\t\t {:.4f}'.format(np.std(hp_cycle)))
except:
print('Error: hp_filter takes only numbers as parameters')
return None
def stoch_simul(self, irf=40):
self.TimeIteration()
FT = self.F.T
QT = self.Q.T
print('\n'+'POLICY AND TRANSITION FUNCTIONS'+'\n')
header = '\t'
for v in self.var_symbols:
header += '\t' + str(v)
print(header)
line = ''
for i in range(self.n):
line += '\t%.4f' % self.steady_values[self.steadys[i]]
print('Constant' + line)
for i in range(self.n):
if (FT[i] != np.zeros((1, self.n))).any():
line = '\t'
for j in range(self.n):
line += '\t%.4f' % FT[i, j]
print(str(self.var_symbols[i])+'(-1)', line)
for i in range(self.p):
line = '\t'
for j in range(self.n):
line += '\t%.4f' % QT[i, j]
print(str(self.varexo_symbols[i]), ' ', line)
# Impulse response functions
if irf > 0:
x = np.zeros((self.n, irf+2))
É› = np.zeros((self.p, irf+2))
É›[:, 1] = 1
for t in range(1, irf+2):
x[:, t] = self.F @ x[:, t-1] + self.Q @ É›[:, t]
for i in range(self.n):
plt.plot(x[i, 1:].T, 'k')
plt.hlines(0, 0, irf, 'r')
plt.title(str(self.var_symbols[i]))
plt.show()
return None
In [30]:
var = 'y c i k h w r z yh R'
varexo = 'e'
param_values = {'alpha':0.33, sy.symbols('beta'):0.99, 'delta':0.025, 'phi':1.75, 'sigma':1, 'rho':0.962159059423142}
model = ('-log(z) + rho*log(z(-1)) + e',
'-k + i + (1-delta)*k(-1)',
'-c^(-sigma) + betta*c(+1)**(-sigma)*(1+r(+1))',
'-y + z*k(-1)^alpha*h^(1-alpha)',
'-r + alpha*y/k(-1) - delta',
'-w + (1-alpha)*y/h',
'-phi/(1-h) + w/c',
'-y + c + i',
'-yh + y/h',
'-R + 1+r')
initval = (1, 0.8, 0.2, 10, 0.33, 2, 0.01, 1, 3, 1.01)
In [31]:
rbc = Dynare(var, varexo, param_values, model, initval)
rbc.system
Out[31]:
In [32]:
rbc.steady()
rbc.stoch_simul(irf=40)
In [33]:
rbc.SimulatedMoments(hp_filter=1600, shocks_stderr=0.007828026718977966)
In [34]:
print('Standard Deviations')
print(hp_cycles.std())