This notebook creates the figures and tables of "Using the Sequence-Space Jacobian to Solve and Estimate Heterogeneous-Agent Models".
You first need to run the notebook models
to get the inputs.
Solution from Reiter method are generated by the notebooks in import_export/autodiff
folder.
reiter_ks.ipynb
reiter_hank1.ipynb
Jacobians with automatic differentiation are generated by the notebooks in import_export/autodiff
folder.
jac_accuracy_ks.ipynb
jac_accuracy_hank.ipynb
# Python packages
import numpy as np
import pandas as pd
from scipy.io import loadmat
import seaborn as sns
import matplotlib.pyplot as plt
import inspect
# Models
from models import krusell_smith as ks
from models import hank_1a
from models import hank_2a
# Auxiliary functions
from toolkit import aux_fn as aux
from toolkit import aux_speed as aux_jac
# Figure fonts
plt.rc('text', usetex=True)
plt.rc('font', family='serif')
plt.rc('text.latex', preamble=r'\usepackage{mathpazo}')
plt.rc('font', size=15)
# Save figures?
savefig = False
# Color palette for optimal black-and-white print
cpalette = sns.color_palette("cubehelix", 5)
sns.palplot(cpalette)
sns.set_palette("cubehelix", 5)
Fun fact: these colors are called regal blue, san felix, muddy waters, plum, pale cornflower blue.
# Load model
ss, G = np.load('import_export/results/hd_ks.npy', allow_pickle=True)
T = 300
Tplot = 51
# Fig 1a
rhos = np.array([0.3, 0.5, 0.7, 0.8, 0.9])
dZ = rhos ** (np.arange(T)[:, np.newaxis])
dK = G['K']['Z'] @ dZ
plt.figure(figsize=(6, 4.5))
for s in range(len(rhos)): plt.plot(dK[:Tplot, s], linewidth=2, label=r'$\rho$=' + str(rhos[s]))
plt.ylabel(r'\% deviation from ss')
plt.xlabel(r'Time $t$')
plt.axhline(y=0, color='#808080', linestyle=':')
plt.legend(framealpha=0)
plt.xlim([0, Tplot - 1])
plt.tight_layout()
if savefig: plt.savefig('import_export/figures/fig1_ks_rho.pdf', format='pdf')
# Compute impulse responses
news = np.array([5, 10, 15, 20, 25])
dZ = np.arange(T)[:, np.newaxis] == news
dK = G['K']['Z'] @ dZ
# Fig 1b
dZ = np.arange(T)[:, np.newaxis] == news
dK = G['K']['Z'] @ dZ
plt.figure(figsize=(6, 4.5))
for s in range(len(news)): plt.plot(dK[:Tplot, s], linewidth=2, label=r'$s$=' + str(news[s]))
plt.ylabel(r'\% deviation from ss')
plt.xlabel(r'Time $t$')
plt.axhline(y=0, color='#808080', linestyle=':')
plt.legend(framealpha=0)
plt.xlim([0, Tplot - 1])
plt.ylim([-.4, 0.9])
plt.tight_layout()
if savefig: plt.savefig('import_export/figures/fig1_ks_news.pdf', format='pdf')
colors = plt.rcParams['axes.prop_cycle'].by_key()['color']
# Load model
T = 300
ss,G = np.load('import_export/results/ks.npy',allow_pickle=True)
J, F = aux_jac.all_JFs(ks.backward_iterate, ss, T, {'r': {'r': 1}, 'w': {'w': 1}})
# Fake news matrix and Jacobian for Krusell-Smith."""
ss = ks.ks_ss()
news=np.array([0, 25, 50, 75, 100])
# fake news col 0
plt.figure(figsize =(4, 3))
plt.plot(F['A']['r'][:, 0], linewidth=2, label=r'$s$=0')
plt.xlabel(r'Time $t$')
plt.xticks(np.arange(0, T+1, 100))
plt.axhline(y=0, color='#808080', linestyle=':')
plt.legend(framealpha=0)
plt.xlim([0,T])
plt.tight_layout()
if savefig: plt.savefig('import_export/figures/fig2_fake0.pdf', format='pdf')
# fake news cols > 0
plt.figure(figsize =(4, 3))
colors = plt.rcParams['axes.prop_cycle'].by_key()['color']
plt.gca().set_prop_cycle(color=colors[1:len(news)])
for s in news[1:]:
plt.plot(F['A']['r'][:, s], linewidth=2, label=r'$s$=' + str(s))
plt.xlabel(r'Time $t$')
plt.xticks(np.arange(0, T+1, 100))
plt.axhline(y=0, color='#808080', linestyle=':')
plt.legend(framealpha=0)
plt.yticks([0,0.1,0.2])
plt.xlim([0,T])
plt.tight_layout()
if savefig: plt.savefig('import_export/figures/fig2_fake_rest.pdf', format='pdf')
# jacobian
plt.figure(figsize =(4, 3))
for s in news:
plt.plot(J['A']['r'][:, s], linewidth=2, label=r'$s$=' + str(s))
plt.xlabel(r'Time $t$')
plt.xticks(np.arange(0, T+1, 100))
plt.yticks(np.arange(0, 13, 4)) # hardcoded!
plt.axhline(y=0, color='#808080', linestyle=':')
plt.xlim([0,T])
plt.legend(framealpha=0)
plt.tight_layout()
if savefig: plt.savefig('import_export/figures/fig2_jac.pdf', format='pdf')
sns.set_palette("cubehelix", 4)
cpalette = sns.color_palette("cubehelix", 4)
# Load series
irf_reiter = pd.read_csv('import_export/autodiff/impulse_responses/irf_reiter_ks.csv', header=None).values
irf_ssj_autodiff = pd.read_csv('import_export/autodiff/impulse_responses/irf_ssj_autodiff_ks.csv', header=None).values
irf_ssj_numdiff = pd.read_csv('import_export/autodiff/impulse_responses/irf_ssj_numdiff_ks.csv', header=None).values
irf_ssj_2sided = pd.read_csv('import_export/autodiff/impulse_responses/irf_ssj_twosidednumdiff_ks.csv', header=None).values
# Compute distance from Reiter
irf_dist_autodiff = np.abs(irf_ssj_autodiff - irf_reiter)
irf_dist_numdiff = np.abs(irf_ssj_numdiff - irf_reiter)
irf_dist_2sided = np.abs(irf_ssj_2sided - irf_reiter)
# Plots
plt.figure(figsize=(6, 4.5))
plt.plot(100 * irf_reiter, label='Reiter', linewidth=2)
plt.plot(100 * irf_ssj_numdiff, label='SSJ with 1-sided num. diff.', linewidth=2, linestyle=':')
plt.plot(100 * irf_ssj_2sided, label='SSJ with 2-sided num. diff.', linewidth=2, linestyle='--')
plt.plot(100 * irf_ssj_autodiff, label='SSJ with auto. diff.', linewidth=2, linestyle='-.')
plt.ylabel(r'\% deviation from ss')
plt.xlabel(r'Time $t$')
plt.xlim([0, len(irf_reiter) - 1])
plt.legend(framealpha=0, loc='upper right')
plt.tight_layout()
if savefig: plt.savefig(f'import_export/figures/fig5_reiter_ks1.pdf', format='pdf')
plt.figure(figsize=(6, 4.5))
plt.plot(100 * irf_dist_numdiff, label='SSJ with 1-sided num. diff.', linewidth=2, linestyle=':', color=cpalette[1])
plt.plot(100 * irf_dist_2sided, label='SSJ with 2-sided num. diff.', linewidth=2, linestyle='--', color=cpalette[2])
plt.plot(100 * irf_dist_autodiff, label='SSJ with auto. diff.', linewidth=2, linestyle='-.', color=cpalette[3])
plt.xlabel(r'Time $t$')
plt.yscale('log')
plt.xlim([0, len(irf_reiter) - 1])
plt.tight_layout()
plt.ylim(1E-12, 1E-1)
plt.legend(framealpha=0, loc='upper right')
if savefig: plt.savefig(f'import_export/figures/fig5_reiter_ks2.pdf', format='pdf')
# Load series
irf_reiter = pd.read_csv('import_export/autodiff/impulse_responses/irf_reiter_hank.csv', header=None).values
irf_ssj_autodiff = pd.read_csv('import_export/autodiff/impulse_responses/irf_ssj_autodiff_hank.csv', header=None).values
irf_ssj_numdiff = pd.read_csv('import_export/autodiff/impulse_responses/irf_ssj_numdiff_hank.csv', header=None).values
irf_ssj_2sided = pd.read_csv('import_export/autodiff/impulse_responses/irf_ssj_twosidednumdiff_hank.csv', header=None).values
# Compute distance from Reiter
irf_dist_autodiff = np.abs(irf_ssj_autodiff - irf_reiter)
irf_dist_numdiff = np.abs(irf_ssj_numdiff - irf_reiter)
irf_dist_2sided = np.abs(irf_ssj_2sided - irf_reiter)
# Plots
plt.figure(figsize=(6, 4.5))
plt.plot(100 * irf_reiter, label='Reiter', linewidth=2)
plt.plot(100 * irf_ssj_numdiff, label='SSJ with 1-sided num. diff.', linewidth=2, linestyle=':')
plt.plot(100 * irf_ssj_2sided, label='SSJ with 2-sided num. diff.', linewidth=2, linestyle='--')
plt.plot(100 * irf_ssj_autodiff, label='SSJ with auto. diff.', linewidth=2, linestyle='-.')
plt.ylabel(r'\% deviation from ss')
plt.xlabel(r'Time $t$')
plt.xlim([0, len(irf_reiter) - 1])
plt.legend(framealpha=0, loc='upper right')
plt.tight_layout()
if savefig: plt.savefig(f'import_export/figures/fig5_reiter_ha1.pdf', format='pdf')
plt.figure(figsize=(6, 4.5))
plt.plot(100 * irf_dist_numdiff, label='SSJ with 1-sided num. diff.', linewidth=2, linestyle=':', color=cpalette[1])
plt.plot(100 * irf_dist_2sided, label='SSJ with 2-sided num. diff.', linewidth=2, linestyle='--', color=cpalette[2])
plt.plot(100 * irf_dist_autodiff, label='SSJ with auto. diff.', linewidth=2, linestyle='-.', color=cpalette[3])
plt.xlabel(r'Time $t$')
plt.yscale('log')
plt.xlim([0, len(irf_reiter) - 1])
plt.tight_layout()
plt.ylim(1E-13, 1E-1)
plt.legend(framealpha=0, loc='upper right')
if savefig: plt.savefig(f'import_export/figures/fig5_reiter_ha2.pdf', format='pdf')
sns.set_palette("cubehelix", 3)
# Load data
T_list, rmse_ks_rho90, rmse_hank_rho90, rmse_hank2_rho90, rmse_ks_rho99, rmse_hank_rho99, rmse_hank2_rho99 = np.load(
'import_export/results/truncation.npy', allow_pickle=True)
# Plot transitory shock
plt.figure(figsize=(6, 4.5))
plt.plot(T_list[:-1], rmse_ks_rho90[:-1], linewidth=2, label='Krusell-Smith')
plt.plot(T_list[:-1], rmse_hank_rho90[:-1], linewidth=2, label='One-asset HANK', linestyle='--')
plt.plot(T_list[:-1], rmse_hank2_rho90[:-1], linewidth=2, label='Two-asset HANK', linestyle='-.')
plt.ylabel(r'RMSE in $dY/Y_{ss}$ in first 100 periods')
plt.ylim(1E-9, 1)
plt.xlim(100, 500)
plt.xlabel(r'Truncation horizon')
plt.yscale('log')
plt.legend(framealpha=0, loc='upper right')
plt.tight_layout()
if savefig: plt.savefig(f'import_export/figures/fig6_truncation1.pdf', format='pdf')
# Plot persistent shock
plt.figure(figsize=(6, 4.5))
plt.plot(T_list[:-1], rmse_ks_rho99[:-1], linewidth=2, label='Krusell-Smith')
plt.plot(T_list[:-1], rmse_hank_rho99[:-1], linewidth=2, label='One-asset HANK', linestyle='--')
plt.plot(T_list[:-1], rmse_hank2_rho99[:-1], linewidth=2, label='Two-asset HANK', linestyle='-.')
plt.ylabel(r'RMSE in $dY/Y_{ss}$ in first 100 periods')
plt.ylim(1E-9, 1)
plt.xlim(100, 500)
plt.xlabel(r'Truncation horizon')
plt.yscale('log')
plt.legend(framealpha=0, loc='upper right')
plt.tight_layout()
if savefig: plt.savefig(f'import_export/figures/fig6_truncation2.pdf', format='pdf')
sns.set_palette("cubehelix", 4)
# Load model jacobian
ss, G = np.load('import_export/results/ks.npy', allow_pickle=True)
sigma = 0.02
rho = 0.9
T = 300
lag = 60
## Compute simulations
T_simul = 200
G['Z'] = {};
G['Z']['Z'] = np.identity(T) # add shock itself to G matrix to compute simulated path
outputs = ['Z', 'Y', 'C', 'K'] # define outputs to compute
np.random.seed(2)
simul = aux.simulate(G, ['Z'], outputs, {'Z': rho}, {'Z': sigma}, T, T_simul)
# Compute impulse response
dZ = rho ** (np.arange(T))
dY, dC, dK = G['Y']['Z'] @ dZ, G['C']['Z'] @ dZ, G['K']['Z'] @ dZ
dX = np.stack([dZ, dY, dC, dK], axis=1)
# Compute covariance
def all_covariances_oneshock(dX, sigma, T):
dft = np.fft.rfftn(dX, s=(2 * T - 2,), axes=(0,))
total = sigma ** 2 * (dft.conjugate()[:, :, np.newaxis] * dft[:, np.newaxis, :])
return np.fft.irfftn(total, s=(2 * T - 2,), axes=(0,))[:T]
Sigma = all_covariances_oneshock(dX, sigma, T)
# get sd of each series and correlation
sd = np.sqrt(np.diag(Sigma[0, ...]))
correl = (Sigma / sd) / (sd[:, np.newaxis])
## Figures
# format
ls = np.arange(-lag, lag + 1)
corrs_l_positive = correl[:lag + 1, 0, :]
corrs_l_negative = correl[lag:0:-1, :, 0]
corrs_combined = np.concatenate([corrs_l_negative, corrs_l_positive])
# Plot simulations
plt.figure(figsize=(6, 4.5))
plt.plot(100 * simul['Z'] / ss['Z'], linewidth=2, label=r'$dZ/Z$')
plt.plot(100 * simul['Y'] / ss['Y'], linewidth=2, label=r'$dY/Y$')
plt.plot(100 * simul['C'] / ss['C'], linewidth=2, label=r'$dC/C$')
plt.plot(100 * simul['K'] / ss['K'], linewidth=2, label=r'$dK/K$')
plt.legend(framealpha=0, loc='upper right')
plt.ylabel(r'\% deviation from ss')
plt.xlim(-0, T_simul)
plt.xlabel(r'Time $t$')
plt.tight_layout()
if savefig: plt.savefig('import_export/figures/fig8_simul.pdf', format='pdf')
# Plot moments
plt.figure(figsize=(6, 4.5))
plt.plot(ls, corrs_combined[:, 0], linewidth=2, label=r'$dZ$')
plt.plot(ls, corrs_combined[:, 1], linewidth=2, label=r'$dY$')
plt.plot(ls, corrs_combined[:, 2], linewidth=2, label=r'$dC$')
plt.plot(ls, corrs_combined[:, 3], linewidth=2, label=r'$dK$')
plt.legend(framealpha=0, loc='upper right')
plt.xlim(-lag, lag)
plt.xlabel(r'Lag $l$')
plt.tight_layout()
if savefig: plt.savefig('import_export/figures/fig8_autocov.pdf', format='pdf')
sns.set_palette("cubehelix", 3)
# Load solution
td_nonlinear, td_linear, td_newss, ss = np.load('import_export/results/hank_2a_transition.npy', allow_pickle=True)
# Plot transitory shock
Tplot = 11
plt.figure(figsize=(6, 4.5))
plt.plot(100 * td_linear[:Tplot], linewidth=2, label='Linear', linestyle='-')
plt.plot(100 * td_nonlinear[:Tplot], linewidth=2, label='Nonlinear', linestyle='--')
plt.axhline(y=0, color='#808080', linestyle=':')
plt.legend(framealpha=0)
plt.ylabel(r'\% deviation from ss')
plt.xlim(0, 10)
plt.xlabel(r'Time $t$')
plt.xticks([0, 5, 10])
plt.tight_layout()
if savefig: plt.savefig('import_export/figures/fig9_nonlinear.pdf', format='pdf')
# Plot permanent shock
Tplot = 200
plt.figure(figsize=(6, 4.5))
plt.plot(100 * (td_newss['Y'][:Tplot] - ss['Y']) / ss['Y'], label=r'$Y$')
plt.plot(100 * (td_newss['N'][:Tplot] - ss['N']) / ss['N'], label=r'$N$', linestyle='--')
plt.plot(100 * (td_newss['K'][:Tplot] - ss['K']) / ss['K'], label=r'$K$', linestyle='-.')
plt.axhline(y=0, color='#808080', linestyle=':')
plt.legend(framealpha=0)
plt.ylabel(r'\% deviation from initial ss')
plt.xlim(-0.3, Tplot)
plt.xlabel(r'Time $t$')
plt.tight_layout()
if savefig: plt.savefig('import_export/figures/fig9_newss.pdf', format='pdf')
sns.set_palette("cubehelix", 5)
colors = plt.rcParams['axes.prop_cycle'].by_key()['color']
# Load series
for s in ['0', '50', '100']:
# Load series
J_fakenews_1S = pd.read_csv('import_export/autodiff/jacobians/J_fakenews_oneside_ks_' + s + '.csv').values
J_direct_1S = pd.read_csv('import_export/autodiff/jacobians/J_direct_oneside_ks_' + s + '.csv').values
J_fakenews_2S = pd.read_csv('import_export/autodiff/jacobians/J_fakenews_twoside_ks_' + s + '.csv').values
J_direct_2S = pd.read_csv('import_export/autodiff/jacobians/J_direct_twoside_ks_' + s + '.csv').values
J_fakenews_AD = pd.read_csv('import_export/autodiff/jacobians/J_fakenews_ad_ks_' + s + '.csv').values
J_direct_AD = pd.read_csv('import_export/autodiff/jacobians/J_direct_ad_ks_' + s + '.csv').values
# Plot and save graph
plt.figure(figsize=(4, 3))
plt.plot(np.abs(J_fakenews_1S - J_direct_AD), linewidth=2)
plt.plot(np.abs(J_direct_1S - J_direct_AD), linewidth=2, linestyle='--')
plt.plot(np.abs(J_fakenews_2S - J_direct_AD), linewidth=2)
plt.plot(np.abs(J_direct_2S - J_direct_AD), linewidth=2, linestyle='--')
plt.plot(np.abs(J_fakenews_AD - J_direct_AD), linewidth=2)
if s == '0':
plt.ylabel(r'vs Direct, auto. diff.')
plt.xlabel(r'Time $t$')
plt.yscale('log')
plt.xlim([0, T - 1])
plt.tight_layout()
if savefig: plt.savefig('import_export/figures/figC1_ks_' + s + '.pdf', format='pdf')
# Now deal with legends
label_list = ['Fake-news, 1-sided', 'Direct, 1-sided', 'Fake-news, 2-sided', 'Direct, 2-sided',
'Fake-news, auto. diff.']
plt.figure(figsize=(11, 1))
plt.plot(np.abs(J_fakenews_1S - J_direct_AD), linewidth=2, label=label_list[0])
plt.plot(np.abs(J_direct_1S - J_direct_AD), linewidth=2, label=label_list[1], linestyle='--')
plt.plot(np.abs(J_fakenews_2S - J_direct_AD), linewidth=2, label=label_list[2])
plt.plot(np.abs(J_direct_2S - J_direct_AD), linewidth=2, label=label_list[3], linestyle='--')
plt.plot(np.abs(J_fakenews_AD - J_direct_AD), linewidth=2, label=label_list[4])
plt.ylim(1, 2)
plt.axis('off')
plt.legend(framealpha=0, loc='center', ncol=3);
if savefig: plt.savefig(f'import_export/figures/figC1_ks_title.pdf', format='pdf');
Need outputs from automatic differentiation in autodiff/jacobians
folder.
for s in ['0', '50', '100']:
# Load series
J_fakenews_1S = pd.read_csv('import_export/autodiff/jacobians/J_fakenews_oneside_hank_' + s + '.csv').values
J_direct_1S = pd.read_csv('import_export/autodiff/jacobians/J_direct_oneside_hank_' + s + '.csv').values
J_fakenews_2S = pd.read_csv('import_export/autodiff/jacobians/J_fakenews_twoside_hank_' + s + '.csv').values
J_direct_2S = pd.read_csv('import_export/autodiff/jacobians/J_direct_twoside_hank_' + s + '.csv').values
J_fakenews_AD = pd.read_csv('import_export/autodiff/jacobians/J_fakenews_ad_hank_' + s + '.csv').values
J_direct_AD = pd.read_csv('import_export/autodiff/jacobians/J_direct_ad_hank_' + s + '.csv').values
# Plot and save graph
plt.figure(figsize=(4, 3))
plt.plot(np.abs(J_fakenews_1S - J_direct_AD), linewidth=2)
plt.plot(np.abs(J_direct_1S - J_direct_AD), linewidth=2, linestyle='--')
plt.plot(np.abs(J_fakenews_2S - J_direct_AD), linewidth=2)
plt.plot(np.abs(J_direct_2S - J_direct_AD), linewidth=2, linestyle='--')
plt.plot(np.abs(J_fakenews_AD - J_direct_AD), linewidth=2)
if s == '0':
plt.ylabel(r'vs Direct, auto. diff.')
plt.xlabel(r'Time $t$')
plt.yscale('log')
plt.xlim([0, T - 1])
plt.tight_layout()
if savefig: plt.savefig('import_export/figures/figC1_hank_' + s + '.pdf', format='pdf')
# Now deal with legends
label_list = ['Fake-news, 1-sided', 'Direct, 1-sided', 'Fake-news, 2-sided', 'Direct, 2-sided',
'Fake-news, auto. diff.']
plt.figure(figsize=(11, 1))
plt.plot(np.abs(J_fakenews_1S - J_direct_AD), linewidth=2, label=label_list[0])
plt.plot(np.abs(J_direct_1S - J_direct_AD), linewidth=2, label=label_list[1], linestyle='--')
plt.plot(np.abs(J_fakenews_2S - J_direct_AD), linewidth=2, label=label_list[2])
plt.plot(np.abs(J_direct_2S - J_direct_AD), linewidth=2, label=label_list[3], linestyle='--')
plt.plot(np.abs(J_fakenews_AD - J_direct_AD), linewidth=2, label=label_list[4])
plt.ylim(1, 2)
plt.axis('off')
plt.legend(framealpha=0, loc='center', ncol=3);
if savefig: plt.savefig(f'import_export/figures/figC1_hank_title.pdf', format='pdf')
sns.set_palette("cubehelix", 7)
# Load data from python
x, diff_irf, irf, irf_dist = np.load('import_export/results/sw.npy', allow_pickle=True)
shocks_list = ['eps_a', 'eps_b', 'eps_g', 'eps_I', 'eps_i', 'eps_p', 'eps_w']
label_list = ['Productivity', 'Monetary', 'Gov. spending', 'Investment', 'Discount factor', 'Price markup',
'Wage markup']
# swap order of discount factor shock and MP shock to be consistent with Herbst-Schorfheide
irf[[4, 1]] = irf[[1, 4]]
irf_dist[[4, 1]] = irf_dist[[1, 4]]
plt.figure(figsize=(6, 4.5))
for k in range(len(shocks_list)): plt.plot(irf[k], label=label_list[k], linewidth=2)
plt.axhline(y=0, color='#808080', linewidth=2, linestyle=':')
plt.ylabel(r'\% deviation from ss')
plt.xlabel(r'Time $t$')
plt.xlim(0, 50)
plt.tight_layout()
if savefig: plt.savefig(f'import_export/figures/figB3_sw1.pdf', format='pdf')
plt.figure(figsize=(6, 4.5))
for k in range(len(shocks_list)): plt.plot(irf_dist[k], label=label_list[k], linewidth=2)
plt.ylabel(r'\% deviation from ss')
plt.xlabel(r'Time $t$')
plt.ylim(1E-15, 1E-3)
plt.xlim(0, 50)
plt.yscale('log')
plt.tight_layout()
if savefig: plt.savefig(f'import_export/figures/figB3_sw2.pdf', format='pdf')
plt.figure(figsize=(11, 1))
for k in range(len(shocks_list)): plt.plot(irf_dist[k], label=label_list[k], linewidth=2)
plt.ylim(1, 2)
plt.axis('off')
plt.legend(framealpha=0, loc='center', ncol=4);
if savefig: plt.savefig(f'import_export/figures/figB3_sw3.pdf', format='pdf')
# Load data from python
sns.set_palette("cubehelix", 3)
x,diff_irf,irf,irf_dist = np.load('import_export/results/hs.npy',allow_pickle=True)
label_list = ['Productivity', 'Monetary','Gov. spending']
shocks_list = ['y_epsz','y_epsR','y_epsg']
plt.figure(figsize =(6, 4.5))
for k in range(len(shocks_list)): plt.plot(irf[k], label = label_list[k], linewidth=2)
plt.axhline(y=0, color='#808080', linewidth=2, linestyle=':')
plt.ylabel(r'\% deviation from ss')
plt.xlabel(r'Time $t$')
plt.xlim(0,50)
plt.tight_layout()
if savefig: plt.savefig(f'import_export/figures/figB3_hs1.pdf', format='pdf');
plt.figure(figsize =(6, 4.5))
for k in range(len(shocks_list)): plt.plot(irf_dist[k], label = label_list[k], linewidth=2)
plt.ylabel(r'\% deviation from ss')
plt.xlabel(r'Time $t$')
#plt.ylim(1E-15, 1E-3)
plt.xlim(0,50)
plt.yscale('log')
plt.tight_layout()
if savefig: plt.savefig(f'import_export/figures/figB3_hs2.pdf', format='pdf');
plt.figure(figsize =(11, 1))
for k in range(len(shocks_list)): plt.plot(irf_dist[k], label = label_list[k], linewidth=2)
plt.ylim(1,2)
plt.axis('off')
plt.legend(framealpha=0,loc='center',ncol=4);
if savefig: plt.savefig(f'import_export/figures/figB3_hs3.pdf', format='pdf');
sns.set_palette("cubehelix", 5)
plt.rc('font', size=20)
# Load data
x, para_avg, para_5, para_95, Nsim, acceptancerate, rmean, kernel_support, kernel_density, para_sim = np.load(
'import_export/results/ks_est.npy', allow_pickle=True)
# Plot
plt_title = ['TFP shock s.d.', 'TFP shock AR-1', 'TFP shock AR-2']
plt.rcParams.update({'figure.max_open_warning': 0})
for i in range(para_avg.shape[0]):
# recursive means
plt.figure(figsize=(.8 * 6, .8 * 4.5))
plt.plot(rmean[i], linewidth=2)
plt.ylabel(r'Parameter value')
plt.xlabel(r'Draw')
plt.title(plt_title[i])
plt.tight_layout()
if savefig: plt.savefig(f'import_export/figures/fig_MH_ks_avg{i}.pdf', format='pdf')
# posterior
plt.figure(figsize=(.8 * 6, .8 * 4.5))
plt.plot(kernel_support[i], kernel_density[i], linewidth=2, label='Posterior density')
plt.axvline(np.mean(para_sim[i]), linewidth=2, color=cpalette[1], linestyle=':', label='Recursive average')
plt.axvline(x[i], linewidth=2, color=cpalette[2], linestyle=':', label='Mode from optimization')
plt.xlabel(r'Parameter value')
plt.ylabel(r'Density')
plt.title(plt_title[i])
plt.xlim(kernel_support[i][0], kernel_support[i][-1])
plt.tight_layout()
if savefig: plt.savefig(f'import_export/figures/fig_MH_ks_density{i}.pdf', format='pdf')
# Create legend as well
plt.figure(figsize=(11, 1))
plt.axvline(linewidth=2, label='Posterior density')
plt.axvline(linewidth=2, color=cpalette[1], linestyle=':', label='Recursive average')
plt.axvline(color=cpalette[2], linestyle=':', label='Mode from optimization')
plt.legend(framealpha=0, loc='upper right')
plt.xlim(1, 2)
plt.axis('off')
plt.legend(framealpha=0, loc='center', ncol=4);
if savefig: plt.savefig(f'import_export/figures/fig_MH_legend.pdf', format='pdf')
# Load data
x, para_avg, para_5, para_95, Nsim, acceptancerate, rmean, kernel_support, kernel_density, para_sim = np.load(
'import_export/results/hank_1a_est.npy', allow_pickle=True)
# Plot
plt_title = ['Monetary shock s.d.', 'Monetary shock AR-1', 'G shock s.d.', 'G shock AR-1', 'P markup shock s.d.',
'P markup shock AR-1']
plt.rcParams.update({'figure.max_open_warning': 0})
for i in range(para_avg.shape[0]):
# recursive means
plt.figure(figsize=(.8 * 6, .8 * 4.5))
plt.plot(rmean[i], linewidth=2)
plt.ylabel(r'Parameter value')
plt.xlabel(r'Draw')
plt.title(plt_title[i])
plt.tight_layout()
if savefig: plt.savefig(f'import_export/figures/fig_MH_hank_1a_avg{i}.pdf', format='pdf')
# posterior
plt.figure(figsize=(.8 * 6, .8 * 4.5))
plt.plot(kernel_support[i], kernel_density[i], linewidth=2, label='Posterior density')
plt.axvline(np.mean(para_sim[i]), linewidth=2, color=cpalette[1], linestyle=':', label='Recursive average')
plt.axvline(x[i], linewidth=2, color=cpalette[2], linestyle=':', label='Mode from optimization')
plt.xlabel(r'Parameter value')
plt.ylabel(r'Density')
plt.title(plt_title[i])
plt.xlim(kernel_support[i][0], kernel_support[i][-1])
plt.tight_layout()
if savefig: plt.savefig(f'import_export/figures/fig_MH_hank_1a_density{i}.pdf', format='pdf')
# Load data
x, para_avg, para_5, para_95, Nsim, acceptancerate, rmean, kernel_support, kernel_density, para_sim = np.load(
'import_export/results/hank_1afull_est.npy', allow_pickle=True)
# Plot
plt_title = ['Monetary shock s.d.', 'Monetary shock AR-1', 'G shock s.d.', 'G shock AR-1', 'P markup shock s.d.',
'P markup shock AR-1', r'$\phi$', r'$\phi_y$', r'$\kappa$']
plt.rcParams.update({'figure.max_open_warning': 0})
for i in range(para_avg.shape[0]):
# recursive means
plt.figure(figsize=(.8 * 6, .8 * 4.5))
plt.plot(rmean[i], linewidth=2)
plt.ylabel(r'Parameter value')
plt.xlabel(r'Draw')
plt.title(plt_title[i])
plt.tight_layout()
if savefig: plt.savefig(f'import_export/figures/fig_MH_hank_1afull_avg{i}.pdf', format='pdf')
# posterior
plt.figure(figsize=(.8 * 6, .8 * 4.5))
plt.plot(kernel_support[i], kernel_density[i], linewidth=2, label='Posterior density')
plt.axvline(np.mean(para_sim[i]), linewidth=2, color=cpalette[1], linestyle=':', label='Recursive average')
plt.axvline(x[i], linewidth=2, color=cpalette[2], linestyle=':', label='Mode from optimization')
plt.xlabel(r'Parameter value')
plt.ylabel(r'Density')
plt.title(plt_title[i])
plt.xlim(kernel_support[i][0], kernel_support[i][-1])
plt.tight_layout()
if savefig: plt.savefig(f'import_export/figures/fig_MH_hank_1afull_density{i}.pdf', format='pdf')
# Load data
x, para_avg, para_5, para_95, Nsim, acceptancerate, rmean, kernel_support, kernel_density, para_sim = np.load(
'import_export/results/hank_2a_est.npy', allow_pickle=True)
# Plot
plt_title = ['TFP shock s.d.', 'TFP shock AR-1', 'G shock s.d.', 'G shock AR-1', r'$\beta$ shock s.d.',
r'$\beta$ shock AR-1',
'Investment shock s.d.', 'Investment shock AR-1', 'Monetary shock s.d.', 'Monetary shock AR-1',
'P markup shock s.d.', 'P markup shock AR-1', 'W markup shock s.d.', 'W markup shock AR-1']
plt.rcParams.update({'figure.max_open_warning': 0})
for i in range(para_avg.shape[0]):
# recursive means
plt.figure(figsize=(.8 * 6, .8 * 4.5))
plt.plot(rmean[i], linewidth=2)
plt.ylabel(r'Parameter value')
plt.xlabel(r'Draw')
plt.title(plt_title[i])
plt.tight_layout()
if savefig: plt.savefig(f'import_export/figures/fig_MH_hank_2a_avg{i}.pdf', format='pdf')
# posterior
plt.figure(figsize=(.8 * 6, .8 * 4.5))
plt.plot(kernel_support[i], kernel_density[i], linewidth=2, label='Posterior density')
plt.axvline(np.mean(para_sim[i]), linewidth=2, color=cpalette[1], linestyle=':', label='Recursive average')
plt.axvline(x[i], linewidth=2, color=cpalette[2], linestyle=':', label='Mode from optimization')
plt.xlabel(r'Parameter value')
plt.ylabel(r'Density')
plt.title(plt_title[i])
plt.xlim(kernel_support[i][0], kernel_support[i][-1])
plt.tight_layout()
if savefig: plt.savefig(f'import_export/figures/fig_MH_hank_2a_density{i}.pdf', format='pdf')
# Load data
x, para_avg, para_5, para_95, Nsim, acceptancerate, rmean, kernel_support, kernel_density, para_sim = np.load(
'import_export/results/hank_2afull_est.npy', allow_pickle=True)
# Plot
plt_title = ['TFP shock s.d.', 'TFP shock AR-1', 'G shock s.d.', 'G shock AR-1', r'$\beta$ shock s.d.',
r'$\beta$ shock AR-1',
'Investment shock s.d.', 'Investment shock AR-1', 'Monetary shock s.d.', 'Monetary shock AR-1',
'P markup shock s.d.', 'P markup shock AR-1', 'W markup shock s.d.', 'W markup shock AR-1',
r'$\phi$', r'$\phi_y$', r'$\kappa_p$', r'$\kappa_w$', r'$\epsilon_I$']
plt.rcParams.update({'figure.max_open_warning': 0})
for i in range(para_avg.shape[0]):
# recursive means
plt.figure(figsize=(.8 * 6, .8 * 4.5))
plt.plot(rmean[i], linewidth=2)
plt.ylabel(r'Parameter value')
plt.xlabel(r'Draw')
plt.title(plt_title[i])
plt.tight_layout()
if savefig: plt.savefig(f'import_export/figures/fig_MH_hank_2afull_avg{i}.pdf', format='pdf')
# posterior
plt.figure(figsize=(.8 * 6, .8 * 4.5))
plt.plot(kernel_support[i], kernel_density[i], linewidth=2, label='Posterior density')
plt.axvline(np.mean(para_sim[i]), linewidth=2, color=cpalette[1], linestyle=':', label='Recursive average')
plt.axvline(x[i], linewidth=2, color=cpalette[2], linestyle=':', label='Mode from optimization')
plt.xlabel(r'Parameter value')
plt.ylabel(r'Density')
plt.title(plt_title[i])
plt.xlim(kernel_support[i][0], kernel_support[i][-1])
plt.tight_layout()
if savefig: plt.savefig(f'import_export/figures/fig_MH_hank_2afull_density{i}.pdf', format='pdf')
# Load data
times_fakenews = np.load('import_export/results/fake_news_speed.npy', allow_pickle=True)[0]
times_ha = times_fakenews['fake'][0]
times_irf = np.load('import_export/results/irf_speed.npy', allow_pickle=True)[0]
times_irf['DAG'][0]['hd_ks'] = times_irf['DAG'][0]['ks']
times_G = np.load('import_export/results/G_speed.npy', allow_pickle=True)[0]
times_G['DAG'][0]['hd_ks'] = times_G['DAG'][0]['ks']
times_est = {}
times_est['ks'] = np.load('import_export/results/ks_speed.npy', allow_pickle=True)[0]
times_est['hd_ks'] = np.load('import_export/results/hd_ks_speed.npy', allow_pickle=True)[0]
times_est['ha'] = np.load('import_export/results/hank_1a_speed.npy', allow_pickle=True)[0]
times_est['hafull'] = np.load('import_export/results/hank_1afull_speed.npy', allow_pickle=True)[0]
times_est['ha2'] = np.load('import_export/results/hank_2a_speed.npy', allow_pickle=True)[0]
times_est['ha2full'] = np.load('import_export/results/hank_2afull_speed.npy', allow_pickle=True)[0]
times_nonlinear = np.load('import_export/results/nonlinear_speed.npy', allow_pickle=True)[0]
# Format tables
t1, t2, t3, t4, t5, t6, t7, t8, t9, t10, t11 = '\n', '', '', '', '\n', '', '', '', '', '', '\n'
tss, tsimul = '\n', ''
para_list = ['Krusell-Smith', 'HD Krusell-Smith', 'one-asset HANK', 'two-asset HANK']
id_list = ['ks', 'hd_ks', 'ha', 'ha2']
for n, var in enumerate(id_list):
t1 += id_list[n] + '\t'
tss += ("{:.2f}").format(times_est[var]['ss']) + '\t'
t2 += ("{:.2f}").format(times_ha[var]) + '\t'
t3 += ("{:.4f}").format(times_irf['DAG'][0][var]) + '\t'
t4 += ("{:.4f}").format(times_G['DAG'][0][var]) + '\t'
t5 += ("{:.4f}").format(times_est[var]['lk']) + '\t'
t6 += ("{:.2f}").format(times_est[var]['est']) + '\t'
t7 += ("{:.0f}").format(times_est[var]['mh']) + '\t'
t11 += ("{:.2f}").format(times_nonlinear[var]) + '\t'
tsimul += ("{:.4f}").format(times_est[var]['sim']) + '\t'
t8 = '\n\t' + '\t' + ("{:.3f}").format(times_est['hafull']['lk']) + '\t' + ("{:.3f}").format(
times_est['ha2full']['lk']) + '\t'
t9 = '\t' + '\t' + ("{:.0f}").format(times_est['hafull']['est']) + '\t' + ("{:.0f}").format(
times_est['ha2full']['est']) + '\t'
t10 = '\t' + '\t' + ("{:.0f}").format(times_est['hafull']['mh']) + '\t' + ("{:.0f}").format(
times_est['ha2full']['mh']) + '\t'
# Print results
print('\nComputing times (in seconds)')
print(t1)
print(tss)
print(t2)
print(t3)
print(t4)
print(tsimul)
print(t5)
print(t6)
print(t7)
print(t8)
print(t9)
print(t10)
print(t11)
# Load data
times = np.load('import_export/results/fake_news_speed.npy', allow_pickle=True)[0]
# Format tables
t1, t2, t3, t4, t5, t6, t7, t8, t9 = '\n', '\n', '', '', '\n', '', '', '', ''
para_list = ['Krusell-Smith', 'HD Krusell-Smith', 'one-asset HANK', 'two-asset HANK']
id_list = ['ks', 'hd_ks', 'ha', 'ha2']
for n, var in enumerate(id_list):
t1 += id_list[n] + '\t'
t2 += ("{:.0f}").format(times['direct'][0][var]) + '\t'
t3 += ("{:.0f}").format(times['direct'][1][var]) + '\t'
t4 += ("{:.0f}").format(times['direct'][2][var]) + '\t'
t5 += ("{:.3f}").format(times['fake'][0][var]) + '\t'
t6 += ("{:.3f}").format(times['fake'][1][var]) + '\t'
t7 += ("{:.3f}").format(times['fake'][2][var]) + '\t'
t8 += ("{:.3f}").format(times['fake'][3][var]) + '\t'
t9 += ("{:.3f}").format(times['fake'][4][var]) + '\t'
# Print results
print('\nComputing times (in seconds)')
print(t1)
print(t2)
print(t3)
print(t4)
print(t5)
print(t6)
print(t7)
print(t8)
print(t9)
# Load data
times = {}
times['ks'] = np.load('import_export/results/ks_speed.npy',allow_pickle=True)[0]
times['ha'] = np.load('import_export/results/hank_1a_speed.npy',allow_pickle=True)[0]
times['hafull'] = np.load('import_export/results/hank_1afull_speed.npy',allow_pickle=True)[0]
times['ha2'] = np.load('import_export/results/hank_2a_speed.npy',allow_pickle=True)[0]
times['ha2full'] = np.load('import_export/results/hank_2afull_speed.npy',allow_pickle=True)[0]
# Format tables
t1, t2, t3, t4, t5, t6, t7, t8, t9, t10, t11 = '\t', '\nLk (ms)\t', '1 (ms)\t', '2 (ms)\t', '3 (ms)\t', '\nOpt (s)\t', 'eval\t', '\nMH (s)\t', 'eval\t', 'rate\t', 'c\t'
para_list = ['Krusell-Smith', 'hank shocks', 'hank shocks + para', 'two-asset hank shocks','two-asset hank shocks + para']
id_list = ['ks','ha','hafull','ha2','ha2full']
for n, var in enumerate(id_list):
t1 += id_list[n] + '\t'
t2 += ("{:.3f}").format(1000*times[var]['lk']) + '\t'
t3 += ("{:.3f}").format(1000*times[var]['lk_step1']) + '\t'
t4 += ("{:.3f}").format(1000*times[var]['lk_step2']) + '\t'
t5 += ("{:.3f}").format(1000*times[var]['lk_step3']) + '\t'
t6 += ("{:.2f}").format(times[var]['est']) + '\t'
t7 += ("{:.0f}").format(times[var]['nfev']) + '\t'
t8 += ("{:.2f}").format(times[var]['mh']) + '\t'
t9 += '100K\t'
t10 += ("{:.3f}").format(times[var]['rate']) + '\t'
t11 += ("{:.2f}").format(times[var]['c']) + '\t'
# Print results
print('\nEstimation times\n')
print(t1)
print(t2)
print(t3)
print(t4)
print(t5)
print(t6)
print(t7)
print(t8)
print(t9)
print(t10)
print(t11)
ss = ks.ks_ss()
sig = inspect.signature(ks.ks_ss)
print('r = ', ss['r'])
print('sigma = ', ss['eis']**(-1))
print('alpha = ', ss['alpha'])
print('delta = ', ss['delta'])
print(sig.parameters['rho'])
print(sig.parameters['sigma'])
print('ne = ', len(ss['e_grid']))
print('nk = ', len(ss['a_grid']))
ss = hank_1a.hank_ss()
sig = inspect.signature(hank_1a.hank_ss)
print('beta = ', round(ss['beta'],3) , ' (r =', ss['r'],')')
print('vphi = ', round(ss['vphi'],3))
print('sigma = ', ss['eis']**(-1))
print('nu = ', ss['frisch']**(-1))
print('bbar = ', ss['a_grid'][0])
print(sig.parameters['rho_s'])
print(sig.parameters['sigma_s'])
print('mu = ', ss['mu'])
print('kappa = ', ss['kappa'])
print('B = ', ss['B'])
print('G = ', round(ss['Y']-ss['C'],3))
print('phi = ', ss['phi'])
print('phi = ', 0)
print('ne = ', len(ss['e_grid']))
print('nk = ', len(ss['a_grid']))
ss = hank_2a.hank_ss(noisy=False)
sig = inspect.signature(hank_2a.hank_ss)
print('beta = ', round(ss['beta'],3) , ' (r =', ss['r'],')')
print('sigma = ', ss['eis']**(-1))
print('chi0 = ', ss['chi0'])
print('chi1 = ', round(ss['chi1'],3), ' (B =', round(ss['B'],3),')')
print('chi2 = ', ss['chi2'])
print('bbar = ', ss['a_grid'][0])
print(sig.parameters['rho_z'])
print(sig.parameters['sigma_z'])
print('----------')
print('vphi = ', round(ss['vphi'],3), '( N=',ss['N'],')')
print('nu = ', ss['frisch']**(-1))
print('muw = ', ss['muw'])
print('kappaw = ', ss['kappaw'])
print('----------')
print('Z = ', round(ss['Z'],3), '( Y=',ss['Y'],')')
print('alpha = ', round(ss['alpha'],3), '( K=',ss['K'],')')
print('mup = ', round(ss['mup'],3), '( p+Bg=',round(ss['A']+ss['B'],3),')')
print('delta = ', round(ss['delta'],3))
print('kappap = ', ss['kappap'])
print('----------')
print('omega = ', ss['omega'])
print('----------')
print('tau = ', round(ss['tax'],3))
print('G = ', ss['G'])
print('Bg = ', ss['Bg'])
print('phi = ', ss['phi'])
print('phiy = ',0)
print('----------')
print('ne = ', len(ss['e_grid']))
print('nb = ', len(ss['b_grid']))
print('na = ', len(ss['a_grid']))
# Load data
times = np.load('import_export/results/G_speed.npy',allow_pickle=True)[0]
# Format tables
t1, t2, t3, t4, t5 = '\n', '\n', '', '', ''
para_list = ['Krusell-Smith', 'one-asset HANK', 'two-asset HANK']
id_list = ['ks','ha','ha2']
for n, var in enumerate(id_list):
t1 += id_list[n] + '\t'
t2 += ("{:.1f}").format(1000*times['DAG'][0][var]) + '\t'
t3 += ("{:.1f}").format(1000*times['DAG'][1][var]) + '\t'
t4 += ("{:.1f}").format(1000*times['DAG'][2][var]) + '\t'
t5 += ("{:.1f}").format(1000*times['DAG'][3][var]) + '\t'
# Print results
print('\nComputing times (in ms)')
print(t1)
print(t2)
print(t3)
print(t4)
print(t5)
# Load data
times = np.load('import_export/results/G_speed.npy',allow_pickle=True)[0]
# Format tables
t1, t2, t3 = '\n', '\n', ''
para_list = ['Krusell-Smith', 'one-asset HANK', 'two-asset HANK']
id_list = ['ks','ha','ha2']
for n, var in enumerate(id_list):
t1 += id_list[n] + '\t'
t2 += ("{:.1f}").format(1000*times['DAG'][0][var]) + '\t'
t3 += ("{:.1f}").format(1000*times['no_DAG'][0][var]) + '\t'
# Print results
print('\nComputing times (in ms)')
print(t1)
print(t2)
print(t3)
# Load data
times = np.load('import_export/results/irf_speed.npy',allow_pickle=True)[0]
# Format tables
t1, t2, t3, t4, t5, t6, t7 = '\n', '\n', '\n', '', '', '', ''
para_list = ['Krusell-Smith', 'one-asset HANK', 'two-asset HANK']
id_list = ['ks','ha','ha2']
for n, var in enumerate(id_list):
t1 += id_list[n] + '\t'
# t2 += ("{:.1f}").format(1000*times['no_DAG'][0][var]) + '\t'
t3 += ("{:.1f}").format(1000*times['DAG'][0][var]) + '\t'
t4 += ("{:.1f}").format(1000*times['DAG'][1][var]) + '\t'
t5 += ("{:.1f}").format(1000*times['DAG'][2][var]) + '\t'
t6 += ("{:.1f}").format(1000*times['DAG'][3][var]) + '\t'
t7 += ("{:.1f}").format(1000*times['DAG'][4][var]) + '\t'
# Print results
print('\nComputing times (in seconds)')
print(t3)
print(t4)
print(t5)
print(t6)
print(t7)
# Load data from python
x, diff_irf, irf, irf_dist = np.load('import_export/results/sw.npy', allow_pickle=True)
# Load estimates from dynare
order_dynare = [0, 7, 1, 8, 2, 9, 3, 10, 4, 11, 5, 12, 14, 6, 13, 15]
x_dynare = loadmat('import_export/dynare/smets_wouters/sw_dynarecode_mode')['xparam1'][order_dynare].T[0]
# Format tables
t1, t2, t3, t4, t5, t6 = 'Params\t', 'SSJ\t', 'Dynare\t', '\nParams\t', 'SSJ\t', 'Dynare\t'
para_list = ['sd_a', 'ar_a', 'sd_b', 'ar_b', 'sd_g', 'ar_g', 'sd_I', 'ar_I', 'sd_i', 'ar_i', 'sd_p', 'ar_p', 'ma_p',
'sd_w', 'ar_w', 'ma_w']
for n, var in enumerate(para_list[:10]):
t1 += para_list[n] + '\t'
t2 += ("{:.3f}").format(x[n]) + '\t'
t3 += ("{:.3f}").format(x_dynare[n]) + '\t'
for n, var in enumerate(para_list[10:]):
t4 += para_list[10:][n] + '\t'
t5 += ("{:.3f}").format(x[10:][n]) + '\t'
t6 += ("{:.3f}").format(x_dynare[10:][n]) + '\t'
# Print results
print('\nMode of the posterior distribution\n')
print(t1)
print(t2)
print(t3)
print(t4)
print(t5)
print(t6)
print(
f'\nMaximum distance in irfs between sequence jacobian and dynare across all shocks and series: {np.max(diff_irf):.3}')
print(
f'Maximum distance in estimation between sequence jacobian and dynare: {100 * np.max(np.abs((x_dynare - x) / x)):.2}%')
# Load data
x, diff_irf, irf_python, irf_dist = np.load('import_export/results/hs.npy', allow_pickle=True)
# Load estiamtes from dynare and marlab
x_dynare = loadmat('import_export/dynare/herbst_schorfheide/hs_dynarecode_mode.mat')['xparam1']
# put in the right order
para_list = ['tau', 'kappa', 'psi1', 'psi2', 'sd_r', 'ar_r', 'sd_g', 'ar_g', 'sd_z', 'ar_z', 'cons_r', 'cons_pi', 'cons_y']
order_dynare = [3, 4, 5, 6, 0, 7, 1, 8, 2, 9, 10, 11, 12]
order_python = [6, 7, 8, 9, 0, 5, 1, 2, 3, 4, 10, 11, 12]
estimate_dynare = x_dynare[order_dynare].T[0]
estimate_python = x[order_python]
# Format tables
t1, t2, t3 = 'Params\t', 'SSJ\t', 'Dynare\t'
for n, var in enumerate(para_list):
t1 += para_list[n] + '\t'
t2 += ("{:.4f}").format(estimate_python[n]) + '\t'
t3 += ("{:.4f}").format(estimate_dynare[n]) + '\t'
# Print results
print('\nMode of the posterior distribution\n')
print(t1)
print(t2)
print(t3)
print(
f'\nMaximum distance in irfs between sequence jacobian and dynare across all shocks and series: {np.max(diff_irf):.3}')
print(
f'Maximum distance in estimation between sequence jacobian and dynare: {100 * np.max(np.abs((estimate_dynare - estimate_python) / estimate_python)):.2}%')
# Load data
x, para_avg, para_5, para_95, Nsim, acceptancerate, rmean, kernel_support, kernel_density, para_sim = np.load(
'import_export/results/ks_est.npy', allow_pickle=True)
# Format tables
t1, t2, t3, t4, t5, t6 = '', '', '', '', '', ''
para_list = ['sd', 'ar1', 'ma1']
for n, var in enumerate(para_list):
t1 += para_list[n] + '\t'
t2 += ("{:.3f}").format(x[n]) + '\t'
t3 += para_list[n] + '\t'
t4 += ("{:.3f}").format(para_avg[n]) + '\t'
t5 += ("{:.3f}").format(para_5[n]) + '\t'
t6 += ("{:.3f}").format(para_95[n]) + '\t'
# Print results
print('\nEstimation results:')
print(t1)
print(t2)
print(f'\nNumber of simulations: {Nsim}')
print(f'Acceptance rate: {acceptancerate:.3}')
print('\nRecursive averages:')
print(t3)
print(t4)
print(t5)
print(t6)
# Load data
x, para_avg, para_5, para_95, Nsim, acceptancerate, rmean, kernel_support, kernel_density, para_sim = np.load(
'import_export/results/hank_1a_est.npy', allow_pickle=True)
# Format tables
t1, t2, t3, t4, t5, t6 = '', '', '', '', '', ''
para_list = ['sd_m', 'ar_m', 'sd_g', 'ar_g', 'sd_p', 'ar_p']
for n, var in enumerate(para_list):
t1 += para_list[n] + '\t'
t2 += ("{:.3f}").format(x[n]) + '\t'
t3 += para_list[n] + '\t'
t4 += ("{:.3f}").format(para_avg[n]) + '\t'
t5 += ("{:.3f}").format(para_5[n]) + '\t'
t6 += ("{:.3f}").format(para_95[n]) + '\t'
# Print results
print('\nEstimation results:')
print(t1)
print(t2)
print(f'\nNumber of simulations: {Nsim}')
print(f'Acceptance rate: {acceptancerate:.3}')
print('\nRecursive averages:')
print(t3)
print(t4)
print(t5)
print(t6)
# Load data
x, para_avg, para_5, para_95, Nsim, acceptancerate, rmean, kernel_support, kernel_density, para_sim = np.load(
'import_export/results/hank_1afull_est.npy', allow_pickle=True)
# Format tables
t1, t2, t3, t4, t5, t6 = '', '', '', '', '', ''
para_list = ['sd_m', 'ar_m', 'sd_g', 'ar_g', 'sd_p', 'ar_p', 'phi', 'phi_y', 'kappa']
for n, var in enumerate(para_list):
t1 += para_list[n] + '\t'
t2 += ("{:.3f}").format(x[n]) + '\t'
t3 += para_list[n] + '\t'
t4 += ("{:.3f}").format(para_avg[n]) + '\t'
t5 += ("{:.3f}").format(para_5[n]) + '\t'
t6 += ("{:.3f}").format(para_95[n]) + '\t'
# Print results
print('\nEstimation results:')
print(t1)
print(t2)
print(f'\nNumber of simulations: {Nsim}')
print(f'Acceptance rate: {acceptancerate:.3}')
print('\nRecursive averages:')
print(t3)
print(t4)
print(t5)
print(t6)
# Load data
x, para_avg, para_5, para_95, Nsim, acceptancerate, rmean, kernel_support, kernel_density, para_sim = np.load(
'import_export/results/hank_2a_est.npy', allow_pickle=True)
# Format tables
t1, t2, t3, t4, t5, t6 = '', '', '', '', '', ''
para_list = ['sd_z', 'ar_z', 'sd_g', 'ar_g', 'sd_beta', 'ar_beta', 'sd_inv', 'ar_inv', 'sd_m', 'ar_m', 'sd_p', 'ar_p',
'sd_w', 'ar_w']
for n, var in enumerate(para_list):
t1 += para_list[n] + '\t'
t2 += ("{:.3f}").format(x[n]) + '\t'
t3 += para_list[n] + '\t'
t4 += ("{:.3f}").format(para_avg[n]) + '\t'
t5 += ("{:.3f}").format(para_5[n]) + '\t'
t6 += ("{:.3f}").format(para_95[n]) + '\t'
# Print results
print('\nEstimation results:')
print(t1)
print(t2)
print(f'\nNumber of simulations: {Nsim}')
print(f'Acceptance rate: {acceptancerate:.3}')
print('\nRecursive averages:')
print(t3)
print(t4)
print(t5)
print(t6)
# Load data
x, para_avg, para_5, para_95, Nsim, acceptancerate, rmean, kernel_support, kernel_density, para_sim = np.load(
'import_export/results/hank_2afull_est.npy', allow_pickle=True)
# Format tables
t1, t2, t3, t4, t5, t6, t1bis, t2bis, t3bis, t4bis, t5bis, t6bis = '', '', '', '', '', '', '', '', '', '', '', ''
para_list = ['sd_z', 'ar_z', 'sd_g', 'ar_g', 'sd_beta', 'ar_beta', 'sd_inv', 'ar_inv', 'sd_m', 'ar_m', 'sd_p', 'ar_p',
'sd_w', 'ar_w', '\nphi', 'phi_y', 'kappa_p', 'kappa_w', 'eps_I']
for n, var in enumerate(para_list[:14]):
t1 += para_list[n] + '\t'
t2 += ("{:.3f}").format(x[n]) + '\t'
t3 += para_list[n] + '\t'
t4 += ("{:.3f}").format(para_avg[n]) + '\t'
t5 += ("{:.3f}").format(para_5[n]) + '\t'
t6 += ("{:.3f}").format(para_95[n]) + '\t'
for n, var in enumerate(para_list[14:]):
t1bis += para_list[14:][n] + '\t'
t2bis += ("{:.3f}").format(x[14:][n]) + '\t'
t3bis += para_list[14:][n] + '\t'
t4bis += ("{:.3f}").format(para_avg[14:][n]) + '\t'
t5bis += ("{:.3f}").format(para_5[14:][n]) + '\t'
t6bis += ("{:.3f}").format(para_95[14:][n]) + '\t'
# Print results
print('\nEstimation results:')
print(t1)
print(t2)
print(t1bis)
print(t2bis)
print(f'\nNumber of simulations: {Nsim}')
print(f'Acceptance rate: {acceptancerate:.3}')
print('\nRecursive averages:')
print(t3)
print(t4)
print(t5)
print(t6)
print(t3bis)
print(t4bis)
print(t5bis)
print(t6bis)