Replication step 1 - Results

This notebook generates results for "Using the Sequence-Space Jacobian to Solve and Estimate Heterogeneous-Agent Models" and saves the data to produce tables and figures.

Full notebook can take 6-10 hours to run. Most of that time is taken up by section 5. All sections are self-contained and can be run in isolation (once imports are loaded, of course).

  1. Herbst-Schorfheide
  2. Smets-Wouters
  3. Krusell-Smith
    3.1 Benchmark
    3.2 High-dimensional
  4. One-asset HANK
    4.1 Only shock processes
    4.2 Shocks and parameters
  5. Two-asset HANK
    5.1 Only shock processes
    5.2 Shocks and parameters
  6. Accuracy checks
    6.1 Truncation of Jacobians
    6.2 Direct method and fake news algorithm
  7. Other timing
    7.1 Heterogeneous agent Jacobians
    7.2 General equilibrium Jacobians
    7.3 Linearized impulse responses
    7.4 Nonlinear impulse responses
In [1]:
# Python packages
import time
import numpy as np
import pandas as pd
import statsmodels.api as sm
from scipy.io import loadmat

# Sequence jacobian package
from toolkit import jacobian as jac
from toolkit.solved_block import solved
from toolkit import estimation as ssj_est
from toolkit import nonlinear

# Models
from models import herbst_schorfheide as hs
from models import smets_wouters as sw
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

# Save output of models?
savedata = True

# MCMC chains
Nsim = 100_000  # number of simulations (200_000 for the paper)
Nburn = 50_000  # number of initial periods thrown away (50_000 for the paper)

1. Herbst-Schorfheide

Calibrate steady state and precompute Jacobian.

In [2]:
# Compute steady state
t = time.time()
ss = hs.calibrate_ss()
elapsed_ss0 = time.time() - t

t = time.time()
ss = hs.calibrate_ss()
elapsed_ss = time.time() - t

# Compute model Jacobian with pre-defined parameters
T = 300
unknowns = ['y', 'pi', 'R']
exogenous = ['epsR', 'g', 'z']
targets = ['yres', 'Rres', 'pires']
outputs = ['gdp', 'infl', 'interest']
block_list = [hs.euler, hs.nkpc, hs.taylor, hs.measurement]
t = time.time()
G = jac.get_G(block_list, exogenous, unknowns, targets, T, ss)
elapsed_G0 = time.time() - t

t = time.time()
G = jac.get_G(block_list, exogenous, unknowns, targets, T, ss)
elapsed_G = time.time() - t

# Show times
print(f'\nTime to compute steady state (burn-in): {1000*round(elapsed_ss0,2)}ms')
print(f'Time to compute model jacobian (burn-in): {1000*round(elapsed_G0,2)}ms')
print(f'\nTime to compute steady state: {1000*round(elapsed_ss,2)}ms')
print(f'Time to compute model jacobian: {1000*round(elapsed_G,2)}ms')
Time to compute steady state (burn-in): 0.0ms
Time to compute model jacobian (burn-in): 410.0ms

Time to compute steady state: 0.0ms
Time to compute model jacobian: 30.0ms

Estimate shock parameters. The data for the estimation was taken from the Matlab replication files available from the book's website.

In [3]:
# Load data
data = pd.read_csv('import_export/data/data_hs.csv', index_col=0, sep=';').values

# Parameters of estimation
shock_series = [('epsR', 0), ('g', 1), ('z', 1)]
params = ['rho_r', 'tau', 'kappa', 'psi1', 'psi2', 'data_r', 'data_pi', 'data_gamma']
x_guess = np.array([0.2, 0.2, 0.8, 0.2, 0.8, 0.8, 2, 0.8, 2, 1, 0.5, 3, 0.5])
priors = [('Invgamma_hs', 0.4, 4), ('Invgamma_hs', 1, 4), ('Uniform', 0, 1), ('Invgamma_hs', 0.5, 4),
          ('Uniform', 0, 1), ('Uniform', 0, 1), ('Gamma', 2, 0.5), ('Uniform', 0, 1), ('Gamma', 1.5, 0.25),
          ('Gamma', 0.5, 0.25), ('Gamma', 0.5, 0.5), ('Gamma', 7, 2), ('Normal', 0.4, 0.2)]
bounds = [(0.1, 0.4), (0.3, 0.99), (0.7, 0.999), (0.01, 0.5), (0.7, 0.99), (0.4, 0.9), (1, 5),
          (0.2, 1), (1., 3), (0.1, 2), (0.2, 0.8), (2.5, 4), (0.04, 0.95)]


# Define function to demean the data using the last 3 parameters
def data_demean_f(x, data):
    data_adj = np.zeros_like(data)
    xdta = x[-data.shape[1]:]
    data_adj[:, 0] = data[:, 0] - xdta[2]
    data_adj[:, 1] = data[:, 1] - xdta[1]
    data_adj[:, 2] = data[:, 2] - xdta[0] - xdta[1] - 4 * xdta[2]
    return data_adj


# Store model information for estimation
jac_info = {'T': T, 'block_list': block_list, 'exogenous': exogenous, 'unknowns': unknowns, 'targets': targets,
            'ss': ss}

# Estimate
t = time.time()
result, x_sd, nfev = aux.estimate(outputs, data, x_guess, shock_series, priors, T, params=params, jac_info=jac_info,
                                  bounds=bounds, sd=False, data_demean_f=data_demean_f, tol=1e-12)
x = result.x
elapsed_est0 = time.time() - t

t = time.time()
result, x_sd, nfev = aux.estimate(outputs, data, x_guess, shock_series, priors, T, params=params, jac_info=jac_info,
                                  bounds=bounds, sd=False, data_demean_f=data_demean_f, tol=1e-12)
x = result.x
elapsed_est = time.time() - t

# Show times
print(f'Time for estimation (burn-in): {round(elapsed_est0, 2)}s')
print(f'Time for estimation: {round(elapsed_est, 2)}s')
Time for estimation (burn-in): 61.77s
Time for estimation: 55.79s

Compute impulse responses and compare with Dynare.

In [4]:
# Construct dict with shock process
dshock = {}
dshock['epsR'] = ss['sigma_r'] * 0 ** np.arange(T)
dshock['g'] = ss['sigma_g'] * ss['rho_g'] ** np.arange(T)
dshock['z'] = ss['sigma_z'] * ss['rho_z'] ** np.arange(T)
Tshow = 50

# Build IRF of endogenous variables in python
irf = {}
for k in G:
    irf[k] = {}
    for l in exogenous:
        if l in G[k]:  irf[k][l] = G[k][l] @ dshock[l]

# Load irf from dynare
res_dynare = loadmat('import_export/dynare/herbst_schorfheide/irf_dynare.mat')

# Compute distance between dynare and python for all variables and all shocks
diff_irf = np.zeros(9)
series_python = [['y', 'g'], ['y', 'epsR'], ['y', 'z'], ['pi', 'g'], ['pi', 'epsR'], ['pi', 'z'], ['R', 'g'],
                 ['R', 'epsR'], ['R', 'z']]
series_dynare = ['y_epsg', 'y_epsR', 'y_epsz', 'pi_epsg', 'pi_epsR', 'pi_epsz', 'R_epsg', 'R_epsR', 'R_epsz']
for i in range(9):
    diff_irf[i] = np.sum(
        np.abs(irf[series_python[i][0]][series_python[i][1]][:10] - res_dynare[series_dynare[i]][:10].T))

# Load irf of output from python
series_python = [['y', 'z'], ['y', 'epsR'], ['y', 'g']]
irf_python = np.zeros((3, Tshow))
for i in range(3):
    irf_python[i] = irf[series_python[i][0]][series_python[i][1]][:Tshow]

# Load irf of output from dynare
series_dynare = ['y_epsz', 'y_epsR', 'y_epsg']
irf_dynare = np.zeros((3, Tshow))
for i in range(3):
    irf_dynare[i] = res_dynare[series_dynare[i]][:Tshow].T

# Compute distance between them
irf_dist = np.abs(irf_python - irf_dynare)
In [5]:
# save data
if savedata: np.save('import_export/results/hs', [x, diff_irf, irf_python, irf_dist])
if savedata: np.save('import_export/results/hs_speed', [elapsed_ss, elapsed_G, elapsed_est])

2. Smets-Wouters

Calibrate steady state and precompute Jacobian.

In [6]:
# Compute steady state
t = time.time()
ss, shocks = sw.calibrate_ss()
elapsed_ss0 = time.time() - t

t = time.time()
ss, shocks = sw.calibrate_ss()
elapsed_ss = time.time() - t

# Compute model jacobian
T = 300
unknowns = ['w', 'wf', 'n', 'nf', 'r', 'rf']
shocks_list = ['eps_a', 'eps_b', 'eps_g', 'eps_I', 'eps_i', 'eps_p', 'eps_w']
targets = ['goods_mkt', 'goods_mkt_f', 'fisher', 'w_res', 'mu_p_f', 'mu_w_f']
blocks = [sw.household, sw.household_f, sw.household_markup, sw.firm, sw.firm_f, sw.nkpc_p, sw.nkpc_w, sw.monetary,
          sw.geneq, sw.geneq_f, sw.measurement]
t = time.time()
G = jac.get_G(blocks, shocks_list, unknowns, targets, T, ss)
elapsed_G = time.time() - t

t = time.time()
G = jac.get_G(blocks, shocks_list, unknowns, targets, T, ss)
elapsed_G0 = time.time() - t

# Show times
print(f'\nTime to compute steady state (burn-in): {1000 * round(elapsed_ss0, 2)}ms')
print(f'Time to compute model jacobian (burn-in): {1000 * round(elapsed_G0, 2)}ms')
print(f'\nTime to compute steady state: {1000 * round(elapsed_ss, 2)}ms')
print(f'Time to compute model jacobian: {1000 * round(elapsed_G, 2)}ms')
Time to compute steady state (burn-in): 0.0ms
Time to compute model jacobian (burn-in): 710.0ms

Time to compute steady state: 0.0ms
Time to compute model jacobian: 1030.0ms

Estimate shock parameters. Data for the estimation was taken from the replication files on the AER website of Smets and Wouters (2007).

In [7]:
# Load data
data = (pd.read_csv('import_export/data/data_sw.csv', index_col=0, sep=';').values)[74:, :]

# Remove intercept from series
data_adj = np.zeros_like(data)
data_adj[:, 0] = data[:, 0] - 0.509  # y
data_adj[:, 1] = data[:, 1] - 0.509  # I
data_adj[:, 2] = data[:, 2] - 0.509  # c
data_adj[:, 3] = data[:, 3] - 0.635  # pi
data_adj[:, 4] = data[:, 4] - 1.427  # i
data_adj[:, 5] = data[:, 5] - 1.743  # n
data_adj[:, 6] = data[:, 6] - 0.509  # w

# Define option for estimation
outputs = ['y_dta', 'I_dta', 'c_dta', 'pi_dta', 'i_dta', 'n_dta', 'w_dta']
shock_series = [('eps_a', 1), ('eps_b', 1), ('eps_g', 1), ('eps_I', 1), ('eps_i', 1), ('eps_p', 2), ('eps_w', 2)]
x_guess = [0.4, 0.8, 0.4, 0.8, 0.4, 0.8, 0.4, 0.8, 0.4, 0.8, 0.4, 0.8, 0.8, 0.4, 0.8, 0.8]
priors = [('Invgamma_hs', 0.1, 2), ('Beta', 0.5, 0.20), ('Invgamma_hs', 0.1, 2), ('Beta', 0.5, 0.20),
          ('Invgamma_hs', 0.1, 2), ('Beta', 0.5, 0.20),
          ('Invgamma_hs', 0.1, 2), ('Beta', 0.5, 0.20), ('Invgamma_hs', 0.1, 2), ('Beta', 0.5, 0.20),
          ('Invgamma_hs', 0.1, 2), ('Beta', 0.5, 0.20), ('Beta', 0.5, 0.20),
          ('Invgamma_hs', 0.1, 2), ('Beta', 0.5, 0.20), ('Beta', 0.5, 0.20)]
bounds = [(0.01, 3), (0.01, 0.9999), (0.025, 5), (0.01, 0.9999), (0.01, 3), (0.01, 0.9999), (0.01, 3), (0.01, 0.9999),
          (0.01, 3), (0.01, 0.9999), (0.01, 3), (0.01, 0.9999), (0.01, 0.9999), (0.01, 3), (0.01, 0.9999), (0.01, 0.9999)]

# Estimate
t = time.time()
result, x_sd, nfev = aux.estimate(outputs, data_adj, x_guess, shock_series, priors, T, G=G, bounds=bounds)
elapsed_est0 = time.time() - t

t = time.time()
result, x_sd, nfev = aux.estimate(outputs, data_adj, x_guess, shock_series, priors, T, G=G, bounds=bounds)
elapsed_est = time.time() - t
x = result.x

# Show times
print(f'Time for estimation (burn-in): {round(elapsed_est0, 2)}s')
print(f'Time for estimation: {round(elapsed_est, 2)}s')
Time for estimation (burn-in): 32.94s
Time for estimation: 32.4s

Compute impulse responses and compare with dynare

In [8]:
# Load irfs from dynare
irf_dynare = loadmat('import_export/dynare/smets_wouters/irf_dynare.mat')

# Define variables
series = ['y', 'c', 'I', 'n', 'pi', 'w', 'i']
series_dynare = [['y_ea', 'y_eb', 'y_eg', 'y_eqs', 'y_em', 'y_epinf', 'y_ew'],
                 ['c_ea', 'c_eb', 'c_eg', 'c_eqs', 'c_em', 'c_epinf', 'c_ew', ],
                 ['inve_ea', 'inve_eb', 'inve_eg', 'inve_eqs', 'inve_em', 'inve_epinf', 'inve_ew'],
                 ['lab_ea', 'lab_eb', 'lab_eg', 'lab_eqs', 'lab_em', 'lab_epinf', 'lab_ew'],
                 ['pinf_ea', 'pinf_eb', 'pinf_eg', 'pinf_eqs', 'pinf_em', 'pinf_epinf', 'pinf_ew'],
                 ['w_ea', 'w_eb', 'w_eg', 'w_eqs', 'w_em', 'w_epinf', 'w_ew'],
                 ['r_ea', 'r_eb', 'r_eg', 'r_eqs', 'r_em', 'r_epinf', 'r_ew']]

# Compute CIR for all series and all shocks
diff_irf = np.zeros((len(series), len(shocks_list)))
for i in range(len(series)):
    for j in range(len(shocks_list)):
        irf1 = (G[series[i]][shocks_list[j]] @ (
                    shocks[shocks_list[j]][0] * aux.arma_irf(np.array(shocks[shocks_list[j]][1]),
                                                             np.array(shocks[shocks_list[j]][2]), T)))[:50]
        irf2 = irf_dynare[series_dynare[i][j]].T
        diff_irf[i, j] = np.max(np.abs(irf2 - irf1))

# Compute irf for output
Tshow = 50
irf = np.zeros((len(shocks_list), Tshow))
for k in range(len(shocks_list)):
    dshock = shocks[shocks_list[k]][0] * aux.arma_irf(np.array(shocks[shocks_list[k]][1]),
                                                      np.array(shocks[shocks_list[k]][2]), T)
    irf[k] = (G['y'][shocks_list[k]] @ dshock)[:Tshow]

# Load irfs of output from dynare
res_dynare = loadmat('import_export/dynare/smets_wouters/irf_dynare.mat')
series_dynare = ['y_ea', 'y_eb', 'y_eg', 'y_eqs', 'y_em', 'y_epinf', 'y_ew']
irf_dynare = np.zeros((len(shocks_list), Tshow))
for k in range(len(series_dynare)):
    irf_dynare[k] = res_dynare[series_dynare[k]][:Tshow].T

# Compute distance between them
irf_dist = np.abs(irf - irf_dynare)
In [9]:
# save data
if savedata: np.save('import_export/results/sw', [x, diff_irf, irf, irf_dist])
if savedata: np.save('import_export/results/sw_speed', [elapsed_ss, elapsed_G, elapsed_est])

3. Krusell-Smith

3.1 Benchmark Krusell-Smith

Calibrate steady state and precompute Jacobian.

In [10]:
times = {}

# Steady state
t = time.time()
ss = ks.ks_ss()
elapsed_ss0 = time.time() - t

# Steady state
t = time.time()
ss = ks.ks_ss()
times['ss'] = time.time() - t

# Compute model jacobian G
T = 300

t = time.time()
G = jac.get_G(block_list=[ks.firm, ks.mkt_clearing, ks.household], exogenous=['Z'], unknowns=['K'],
              targets=['asset_mkt'], T=T, ss=ss)
elapsed_G0 = time.time() - t

t = time.time()
G = jac.get_G(block_list=[ks.firm, ks.mkt_clearing, ks.household], exogenous=['Z'], unknowns=['K'],
              targets=['asset_mkt'], T=T, ss=ss)
times['G'] = time.time() - t

# Show times
t1, t2 = times['ss'], times['G']
print(f'\nTime to compute steady state (burn-in): {round(elapsed_ss0, 2)}s')
print(f'Time to compute model jacobian (burn-in): {round(elapsed_G0, 2)}s')

print(f'\nTime to compute steady state: {round(t1, 2)}s')
print(f'Time to compute model jacobian: {round(t2, 2)}s')
Time to compute steady state (burn-in): 1.82s
Time to compute model jacobian (burn-in): 0.48s

Time to compute steady state: 0.41s
Time to compute model jacobian: 0.11s

Bayesian estimation

In [11]:
# Define series to estimate
series = ['y']
data = aux.get_normalized_data(ss, 'import_export/data/data_bayes.csv', series)

# Define option for estimation
shock_series = [('Z', 2)]  # specifies shock to hit economy and whether it's an AR-1 or ARMA 1-1
x_guess = [0.4, 0.5, 0.4]
priors = [('Invgamma', 0.4, 4), ('Beta', 0.5, 0.2), ('Beta', 0.5, 0.2)]
bounds = [(1e-2, 4), (1e-2, 0.98), (1e-2, 0.98)]
outputs = ['Y']

# Estimate
t = time.time()
result, x_sd, times['nfev'] = aux.estimate(outputs, data, x_guess, shock_series, priors, T, G=G, bounds=bounds)
x = result.x
elapsed_est0 = time.time() - t

t = time.time()
result, x_sd, times['nfev'] = aux.estimate(outputs, data, x_guess, shock_series, priors, T, G=G, bounds=bounds)
x = result.x
times['est'] = time.time() - t

# Show times
t1 = times['est']
print(f'Time for estimation (burn-in): {round(t1, 2)}s')
print(f'Time for estimation: {round(t1, 2)}s')
Time for estimation (burn-in): 0.06s
Time for estimation: 0.06s

Simulation at the posterior mode

In [12]:
T_simul = 100_000
t = time.time()
sim = aux.simulate(G, ['Z'], outputs, {'Z': x[0]}, {'Z': x[1]}, T, T_simul)  # Note, just simulate AR(1) for now
times_sim0 = time.time() - t

t = time.time()
sim = aux.simulate(G, ['Z'], outputs, {'Z': x[0]}, {'Z': x[1]}, T, T_simul)  # Note, just simulate AR(1) for now
times['sim'] = time.time() - t

# Show times
t1 = times['sim']
print(f'Time for simulation (burn-in): {round(1000 * times_sim0, 3)}ms')
print(f'Time for simulation: {round(1000 * t1, 3)}ms')
Time for simulation (burn-in): 448.376ms
Time for simulation: 4.987ms

Get run time for likelihood evaluation

In [13]:
T_irf = T - 20
n_se,n_sh = len(outputs), len(shock_series)
meas_error = np.zeros(n_se)
sigmas, arcoefs, macoefs = aux.get_shocks_arma(x, shock_series)

# run them all once (burn-in)
As = aux.step1_est(G, arcoefs, macoefs, shock_series, outputs, T, T_irf, n_se, n_sh)
Sigma = ssj_est.all_covariances(As, sigmas)
llh = aux.step3_est(Sigma, data, sigma_o=meas_error)

# Compute run time
res = %timeit -o -q As = aux.step1_est(G, arcoefs, macoefs, shock_series, outputs, T, T_irf, n_se, n_sh)
times['lk_step1'] = np.mean(res.timings)
res = %timeit -o -q Sigma = ssj_est.all_covariances(As, sigmas)
times['lk_step2'] = np.mean(res.timings)
res = %timeit -o -q llh = aux.step3_est(Sigma, data, sigma_o=meas_error)
times['lk_step3'] = np.mean(res.timings)
times['lk'] = times['lk_step1'] + times['lk_step2'] + times['lk_step3']

Metropolis Hastings

In [14]:
# Parameters
c = 2.5             # scale for the variance
x0 = x.copy()


# Define function that returns log-likelihood + log-prior
def lik_f(x):
    return aux.loglik_f(x, data, outputs, shock_series, priors, T, G)


# Get mode and variance covariance matrix of the log likelihood at the mode
mode = x                                              # mode of the parameters
Sigma = - np.linalg.inv(aux.hessian(lik_f, mode)[0])  # inverse variance-covariance matrix

# MH algorithm
t = time.time()
para_sim, para_avg, para_5, para_95, logposterior, acceptancerate = aux.metropolis_hastings(lik_f, mode, Sigma, bounds,
                                                                                            Nsim, Nburn, c)
times_mh0 = time.time() - t

t = time.time()
para_sim, para_avg, para_5, para_95, logposterior, acceptancerate = aux.metropolis_hastings(lik_f, mode, Sigma, bounds,
                                                                                            Nsim, Nburn, c)
times['mh'] = time.time() - t
times['c'], times['rate'] = c, acceptancerate

# Compute recursive mean
rmean = np.zeros((para_avg.shape[0], para_sim.shape[1]))
for i in range(para_sim.shape[1]):
    rmean[:, i] = np.mean(para_sim[:, :i + 1], 1)

# Compute kernel density
kde = sm.nonparametric.KDEUnivariate(para_sim[0])
kde.fit()
kernel_support = np.zeros((para_avg.shape[0], kde.support.shape[0]))
kernel_density = np.zeros((para_avg.shape[0], kde.support.shape[0]))
for i in range(para_avg.shape[0]):
    kde = sm.nonparametric.KDEUnivariate(para_sim[i])
    kde.fit()
    kernel_support[i] = kde.support
    kernel_density[i] = kde.density

# Show times
t1 = times['mh']
print(f'Time for Metropolis-Hastings algorithm (burn-in): {round(times_mh0, 2)}s')
print(f'Time for Metropolis-Hastings algorithm: {round(t1, 2)}s')
print(f'Acceptance rate: {round(acceptancerate, 2)}')
Time for Metropolis-Hastings algorithm (burn-in): 0.77s
Time for Metropolis-Hastings algorithm: 0.71s
Acceptance rate: 0.26
In [15]:
# save data
if savedata: np.save('import_export/results/ks', [ss, G])
if savedata: np.save('import_export/results/ks_est',
                     [x, para_avg, para_5, para_95, Nsim, acceptancerate, rmean, kernel_support, kernel_density, para_sim])
if savedata: np.save('import_export/results/ks_speed', [times])

3.2 High-dimensional Krusell-Smith

Calibrate steady state and precompute Jacobian.

In [16]:
times = {}

# Steady state burn-in
t = time.time()
ss = ks.ks_ss(nS=50, nA=5000)
elapsed_ss0 = time.time() - t

# Steady state
t = time.time()
ss = ks.ks_ss(nS=50, nA=5000)
times['ss'] = time.time() - t

# Compute model jacobian G
T = 300

t = time.time()
G = jac.get_G(block_list=[ks.firm, ks.mkt_clearing, ks.household], exogenous=['Z'], unknowns=['K'],
              targets=['asset_mkt'], T=T, ss=ss)
elapsed_G0 = time.time() - t

t = time.time()
G = jac.get_G(block_list=[ks.firm, ks.mkt_clearing, ks.household], exogenous=['Z'], unknowns=['K'],
              targets=['asset_mkt'], T=T, ss=ss)
times['G'] = time.time() - t

# Show times
t1, t2 = times['ss'], times['G']
print(f'\nTime to compute steady state (burn-in): {round(elapsed_ss0, 2)}s')
print(f'Time to compute model jacobian (burn-in): {round(elapsed_G0, 2)}s')

print(f'\nTime to compute steady state: {round(t1, 2)}s')
print(f'Time to compute model jacobian: {round(t2, 2)}s')
Time to compute steady state (burn-in): 48.51s
Time to compute model jacobian (burn-in): 11.28s

Time to compute steady state: 48.47s
Time to compute model jacobian: 10.47s

Bayesian estimation

In [17]:
# Define series to estimate
series = ['y']
data = aux.get_normalized_data(ss, 'import_export/data/data_bayes.csv', series)

# Define option for estimation
shock_series = [('Z', 2)]  # specifies shock to hit economy and whether it's an AR-1 or ARMA 1-1
x_guess = [0.4, 0.5, 0.4]
priors = [('Invgamma', 0.4, 4), ('Beta', 0.5, 0.2), ('Beta', 0.5, 0.2)]
bounds = [(1e-2, 4), (1e-2, 0.98), (1e-2, 0.98)]
outputs = ['Y']

# Estimate
t = time.time()
result, x_sd, times['nfev'] = aux.estimate(outputs, data, x_guess, shock_series, priors, T, G=G, bounds=bounds)
elapsed_est0 = time.time() - t

t = time.time()
result, x_sd, times['nfev'] = aux.estimate(outputs, data, x_guess, shock_series, priors, T, G=G, bounds=bounds)
x = result.x
times['est'] = time.time() - t

# Show times
t1 = times['est']
print(f'Time for estimation (burn-in): {round(elapsed_est0, 2)}s')
print(f'Time for estimation: {round(t1, 2)}s')
Time for estimation (burn-in): 0.07s
Time for estimation: 0.07s

Simulation at the posterior mode

In [18]:
T_simul = 100_000
t = time.time()
sim = aux.simulate(G, ['Z'], outputs, {'Z': x[0]}, {'Z': x[1]}, T, T_simul)  # note, just simulate AR(1) for now
elapsed_sim0 = time.time() - t

t = time.time()
sim = aux.simulate(G, ['Z'], outputs, {'Z': x[0]}, {'Z': x[1]}, T, T_simul)  # note, just simulate AR(1) for now
times['sim'] = time.time() - t

# Show times
t1 = times['sim']
print(f'Time for simulation: {round(1000 * elapsed_sim0, 3)}ms')
print(f'Time for simulation (burn-in): {round(1000 * t1, 3)}ms')
Time for simulation: 4.987ms
Time for simulation (burn-in): 4.986ms

Get run time for likelihood evaluation

In [19]:
T_irf = T - 20
n_se,n_sh = len(outputs), len(shock_series)
meas_error = np.zeros(n_se)
sigmas, arcoefs, macoefs = aux.get_shocks_arma(x, shock_series)

# run them all once (burn-in)
As = aux.step1_est(G, arcoefs, macoefs, shock_series, outputs, T, T_irf, n_se, n_sh)
Sigma = ssj_est.all_covariances(As, sigmas)
llh = aux.step3_est(Sigma, data, sigma_o=meas_error)

# Compute run time
res = %timeit -o -q As = aux.step1_est(G, arcoefs, macoefs, shock_series, outputs, T, T_irf, n_se, n_sh)
times['lk_step1'] = np.mean(res.timings)
res = %timeit -o -q Sigma = ssj_est.all_covariances(As, sigmas)
times['lk_step2'] = np.mean(res.timings)
res = %timeit -o -q llh = aux.step3_est(Sigma, data, sigma_o=meas_error)
times['lk_step3'] = np.mean(res.timings)
times['lk'] = times['lk_step1'] + times['lk_step2'] + times['lk_step3']

Metropolis Hastings

In [20]:
# Parameters
c = 2.5             # scale for the variance
x0 = x.copy()

# Define function that returns log-likelihood + log-prior
def lik_f(x):
    return aux.loglik_f(x, data, outputs, shock_series, priors, T, G)

# Get mode and variance covariance matrix of the log likelihood at the mode
mode = x                                              # mode of the parameters
Sigma = - np.linalg.inv(aux.hessian(lik_f, mode)[0])  # inverse variance-covariance matrix

# MH algorithm
t = time.time()
para_sim, para_avg, para_5, para_95, logposterior, acceptancerate = aux.metropolis_hastings(lik_f, mode, Sigma, bounds, Nsim, Nburn, c)
elapsed_mh0 = time.time() - t

t = time.time()
para_sim, para_avg, para_5, para_95, logposterior, acceptancerate = aux.metropolis_hastings(lik_f, mode, Sigma, bounds, Nsim, Nburn, c)
times['mh'] = time.time() - t
times['c'], times['rate'] = c, acceptancerate

# Compute recursive mean
rmean = np.zeros((para_avg.shape[0],para_sim.shape[1]))
for i in range(para_sim.shape[1]):
    rmean[:,i] = np.mean(para_sim[:,:i+1],1)

# Compute kernel density
kde = sm.nonparametric.KDEUnivariate(para_sim[0])
kde.fit()
kernel_support = np.zeros((para_avg.shape[0],kde.support.shape[0]))
kernel_density = np.zeros((para_avg.shape[0],kde.support.shape[0]))
for i in range(para_avg.shape[0]):
    kde = sm.nonparametric.KDEUnivariate(para_sim[i])
    kde.fit()
    kernel_support[i] = kde.support
    kernel_density[i] = kde.density
    
# Show times
t1 = times['mh']
print(f'Time for Metropolis-Hastings algorithm: {round(elapsed_mh0,2)}s')
print(f'Time for Metropolis-Hastings algorithm: {round(t1,2)}s')
print(f'Acceptance rate: {round(acceptancerate,2)}')
C:\Users\Bence\Home\Projects\ssjreplication\main\aux_fn.py:164: RuntimeWarning: invalid value encountered in log
  sum_log_priors += (alpha - 1) * np.log(x[n]) + (beta - 1) * np.log(1 - x[n])
[ 0.18866681  0.93410722 -0.02575396]
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-20-3284d1136d38> in <module>
     13 # MH algorithm
     14 t = time.time()
---> 15 para_sim, para_avg, para_5, para_95, logposterior, acceptancerate = aux.metropolis_hastings(lik_f, mode, Sigma, bounds, Nsim, Nburn, c)
     16 elapsed_mh0 = time.time() - t
     17 

~\Home\Projects\ssjreplication\main\aux_fn.py in metropolis_hastings(lik_f, mode, Sigma, bounds, Nsim, Nburn, c)
    339     para_sim = np.zeros((mode.shape[0], Nsim))
    340     para_sim[:, 0] = x
--> 341     obj = lik_f(x)
    342     logposterior = obj * np.ones(Nsim)
    343     accept = 0

<ipython-input-20-3284d1136d38> in lik_f(x)
      5 # Define function that returns log-likelihood + log-prior
      6 def lik_f(x):
----> 7     return aux.loglik_f(x, data, outputs, shock_series, priors, T, G)
      8 
      9 # Get mode and variance covariance matrix of the log likelihood at the mode

~\Home\Projects\ssjreplication\main\aux_fn.py in loglik_f(x, data, outputs, shock_series, priors_list, T, G)
    231 
    232     # compute the posterior by adding the log prior
--> 233     log_posterior = llh + log_priors(x, priors_list)
    234 
    235     return log_posterior

~\Home\Projects\ssjreplication\main\aux_fn.py in log_priors(x, priors_list)
    168     if np.isinf(sum_log_priors) or np.isnan(sum_log_priors):
    169         print(x)
--> 170         raise ValueError('Need tighter bounds to prevent prior value = 0')
    171     return sum_log_priors
    172 

ValueError: Need tighter bounds to prevent prior value = 0
In [21]:
# only need to save speed data
if savedata: np.save('import_export/results/hd_ks', [ss, G])
if savedata: np.save('import_export/results/hd_ks_est',
                     [x, para_avg, para_5, para_95, Nsim, acceptancerate, rmean, kernel_support, kernel_density, para_sim])
if savedata: np.save('import_export/results/hd_ks_speed', [times])

4. One-asset HANK

4.1 One-asset HANK (only shock processes)

Calibrate steady state and precompute Jacobian.

In [22]:
times = {}

# Steady state burn-in
t = time.time()
ss = hank_1a.hank_ss()
elapsed_ss0 = time.time() - t

# Steady state
t = time.time()
ss = hank_1a.hank_ss()
times['ss'] = time.time() - t

# setup
T = 300
exogenous = ['G', 'markup', 'rstar']
unknowns = ['pi', 'w', 'Y']
targets = ['nkpc_res', 'asset_mkt', 'labor_mkt']

# Compute model jacobian G
block_list = [hank_1a.household_trans, hank_1a.firm, hank_1a.monetary, hank_1a.fiscal, hank_1a.nkpc, hank_1a.mkt_clearing]
t = time.time()
G = jac.get_G(block_list, exogenous, unknowns, targets, T, ss, save=True)
elapsed_G0 = time.time() - t

t = time.time()
G = jac.get_G(block_list, exogenous, unknowns, targets, T, ss, save=True)
times['G'] = time.time() - t

# Show times
t1, t2 = times['ss'], times['G']
print(f'\nTime to compute steady state (burn-in): {round(elapsed_ss0, 2)}s')
print(f'Time to compute model jacobian (burn-in): {round(elapsed_G0, 2)}s')

print(f'\nTime to compute steady state: {round(t1, 2)}s')
print(f'Time to compute model jacobian: {round(t2, 2)}s')
Time to compute steady state (burn-in): 2.49s
Time to compute model jacobian (burn-in): 0.84s

Time to compute steady state: 2.06s
Time to compute model jacobian: 0.71s

Bayesian estimation

In [23]:
# Define series to estimate
series = ['y', 'pi', 'i']
data = aux.get_normalized_data(ss,'import_export/data/data_bayes.csv', series)

# Define option for estimation
shock_series = [('rstar', 1), ('G', 1), ('markup', 1)]
x_guess = [0.4, 0.5] * len(shock_series) 
priors = [('Invgamma', 0.4, 4), ('Beta', 0.5, 0.2)] * len(shock_series)
bounds = [(0.01, 4), (0.01, 0.98)] * len(shock_series)
outputs = ['Y', 'pi', 'i']

# Estimate
t = time.time()
result, x_sd, times['nfev'] = aux.estimate(outputs, data, x_guess, shock_series, priors, T, G=G, bounds=bounds)
x = result.x
times_est0 = time.time() - t

t = time.time()
result, x_sd, times['nfev'] = aux.estimate(outputs, data, x_guess, shock_series, priors, T, G=G, bounds=bounds)
x = result.x
times['est'] = time.time() - t

# Show times
t1 = times['est']
print(f'Time for estimation (burn-in): {round(t1,2)}s')
print(f'Time for estimation: {round(t1,2)}s')
Time for estimation (burn-in): 0.67s
Time for estimation: 0.67s

Simulation at the posterior mode

In [24]:
T_simul = 100_000
t = time.time()
sim = aux.simulate(G, ['rstar', 'G', 'markup'], outputs,
                   {'rstar': x[0], 'G': x[2], 'markup': x[4]},
                   {'rstar': x[1], 'G': x[3], 'markup': x[5]}, T, T_simul)  # note, just simulate AR(1) for now
elapsed_sim0 = time.time() - t

t = time.time()
sim = aux.simulate(G, ['rstar', 'G', 'markup'], outputs,
                   {'rstar': x[0], 'G': x[2], 'markup': x[4]},
                   {'rstar': x[1], 'G': x[3], 'markup': x[5]}, T, T_simul)  # note, just simulate AR(1) for now
times['sim'] = time.time() - t

# Show times
t1 = times['sim']
print(f'Time for simulation (burn-in): {round(1000 * elapsed_sim0, 3)}ms')
print(f'Time for simulation: {round(1000 * t1, 3)}ms')
Time for simulation (burn-in): 22.987ms
Time for simulation: 31.91ms

Get run time for likelihood evaluation

In [25]:
T_irf = T - 20
n_se,n_sh = len(outputs), len(shock_series)
meas_error = np.zeros(n_se)
sigmas, arcoefs, macoefs = aux.get_shocks_arma(x, shock_series)

# run them all once (burn-in)
As = aux.step1_est(G, arcoefs, macoefs, shock_series, outputs, T, T_irf, n_se, n_sh)
Sigma = ssj_est.all_covariances(As, sigmas)
llh = aux.step3_est(Sigma, data, sigma_o=meas_error)

# Compute run time
res = %timeit -o -q As = aux.step1_est(G, arcoefs, macoefs, shock_series, outputs, T, T_irf, n_se, n_sh)
times['lk_step1'] = np.mean(res.timings)
res = %timeit -o -q Sigma = ssj_est.all_covariances(As, sigmas)
times['lk_step2'] = np.mean(res.timings)
res = %timeit -o -q llh = aux.step3_est(Sigma, data, sigma_o=meas_error)
times['lk_step3'] = np.mean(res.timings)
times['lk'] = times['lk_step1'] + times['lk_step2'] + times['lk_step3']

Metropolis Hastings

In [26]:
# Parameters
c = 1.1             # scale for the variance
x0 = x.copy()


# Define function that returns log-likelihood + log-prior
def lik_f(x):
    return aux.loglik_f(x, data, outputs, shock_series, priors, T, G)


# Get mode and variance covariance matrix of the log likelihood at the mode
mode = x                                              # mode of the parameters
Sigma = - np.linalg.inv(aux.hessian(lik_f, mode)[0])  # inverse variance-covariance matrix

# MH algorithm
t = time.time()
para_sim, para_avg, para_5, para_95, logposterior, acceptancerate = aux.metropolis_hastings(lik_f, mode, Sigma, bounds,
                                                                                            Nsim, Nburn, c)
elapsed_mh0 = time.time() - t

t = time.time()
para_sim, para_avg, para_5, para_95, logposterior, acceptancerate = aux.metropolis_hastings(lik_f, mode, Sigma, bounds,
                                                                                            Nsim, Nburn, c)
times['mh'] = time.time() - t
times['c'], times['rate'] = c, acceptancerate

# Compute recursive mean
rmean = np.zeros((para_avg.shape[0], para_sim.shape[1]))
for i in range(para_sim.shape[1]):
    rmean[:, i] = np.mean(para_sim[:, :i + 1], 1)

# Compute kernel density
kde = sm.nonparametric.KDEUnivariate(para_sim[0])
kde.fit()
kernel_support = np.zeros((para_avg.shape[0], kde.support.shape[0]))
kernel_density = np.zeros((para_avg.shape[0], kde.support.shape[0]))
for i in range(para_avg.shape[0]):
    kde = sm.nonparametric.KDEUnivariate(para_sim[i])
    kde.fit()
    kernel_support[i] = kde.support
    kernel_density[i] = kde.density

# Show times
t1 = times['mh']
print(f'Time for Metropolis-Hastings algorithm (burn-in): {round(elapsed_mh0, 2)}s')
print(f'Time for Metropolis-Hastings algorithm: {round(t1, 2)}s')
print(f'Acceptance rate: {round(acceptancerate, 2)}')
Time for Metropolis-Hastings algorithm (burn-in): 2.76s
Time for Metropolis-Hastings algorithm: 2.71s
Acceptance rate: 0.27
In [27]:
# save data
if savedata: np.save('import_export/results/hank_1a_est',
                     [x, para_avg, para_5, para_95, Nsim, acceptancerate, rmean, kernel_support, kernel_density, para_sim])
if savedata: np.save('import_export/results/hank_1a_speed', [times])

4.2. One-asset HANK (shocks + parameters)

Bayesian estimation

In [28]:
times = {}

# Steady state
t = time.time()
ss = hank_1a.hank_ss()
times['ss'] = time.time() - t

# setup
T = 300
exogenous = ['G', 'markup', 'rstar']
unknowns = ['pi', 'w', 'Y']
targets = ['nkpc_res', 'asset_mkt', 'labor_mkt']

# Define series to estimate
series = ['y', 'pi', 'i']
outputs = ['Y', 'pi', 'i']
data = aux.get_normalized_data(ss, 'import_export/data/data_bayes.csv', series)

# Parameters of estimation
shock_series = [('rstar', 1), ('G', 1), ('markup', 1)]
params = ['phi', 'phi_y', 'kappa']
x_guess = [0.4, 0.5] * len(shock_series) + [ss[k] for k in params]
bounds = [(0.01, 4), (0.01, 0.98)] * len(shock_series) + [(1, 3), (1e-3, 1), (1e-4, 0.5)]
priors = [('Invgamma', 0.4, 4), ('Beta', 0.5, 0.2)] * len(shock_series) + [('Gamma', 1.5, 0.25), ('Gamma', 0.5, 0.25),
                                                                           ('Gamma', 0.1, 0.05)]

# Compute jacobian of household block
J_ha = jac.get_G([hank_1a.household_trans], ['r', 'w', 'Div', 'Tax'], [], [], T, ss, save=True)
jac_info = {'T': T, 'J_ha': J_ha,
            'block_list': [J_ha, hank_1a.firm, hank_1a.monetary, hank_1a.fiscal, hank_1a.nkpc, hank_1a.mkt_clearing],
            'exogenous': exogenous, 'unknowns': unknowns, 'targets': targets, 'ss': ss}

# Estimate
t = time.time()
result, x_sd, times['nfev'] = aux.estimate(outputs, data, x_guess, shock_series, priors, T, params=params,
                                           jac_info=jac_info, bounds=bounds, sd=False)
x = result.x
times['est'] = time.time() - t

# Show times
t1, t2 = times['ss'], times['est']
print(f'\nTime to compute steady state: {round(t1, 2)}s')
print(f'Time for estimation: {round(t2, 2)}s')
Time to compute steady state: 2.09s
Time for estimation: 13.35s

Get run time for likelihood evaluation

In [29]:
n_se,n_sh = len(outputs), len(shock_series)
meas_error = np.zeros(n_se)
sigmas, arcoefs, macoefs = aux.get_shocks_arma(x, shock_series)

# run them all once
As = aux.step1_estfull(x,shock_series,outputs,T,params,jac_info)
Sigma = ssj_est.all_covariances(As, sigmas)
llh = aux.step3_est(Sigma, data, sigma_o=meas_error)

# Compute run time
res = %timeit -o -q As = aux.step1_estfull(x,shock_series,outputs,T,params,jac_info)
times['lk_step1'] = np.mean(res.timings)
res = %timeit -o -q Sigma = ssj_est.all_covariances(As, sigmas)
times['lk_step2'] = np.mean(res.timings)
res = %timeit -o -q llh = aux.step3_est(Sigma, data, sigma_o=meas_error)
times['lk_step3'] = np.mean(res.timings)
times['lk'] = times['lk_step1'] + times['lk_step2'] + times['lk_step3']

Metropolis Hastings

In [30]:
# Parameters
c = 0.65            # scale for the variance
x0 = x.copy()


# Define function that returns log-likelihood + log-prior
def lik_f(x):
    n_params = len(params)
    x_params = x[-n_params:]

    # write new parameters into ss
    ss_here = jac_info['ss'].copy()
    ss_here.update({param: x_params[j] for j, param in enumerate(params)})

    # Compute model jacobian G
    G = jac.get_G(jac_info['block_list'], exogenous=jac_info['exogenous'], unknowns=jac_info['unknowns'],
                  targets=jac_info['targets'], T=jac_info['T'], ss=ss_here, outputs=outputs)

    # return likelihood
    return aux.loglik_f(x, data, outputs, shock_series, priors, T, G)


# Get mode and variance covariance matrix of the log likelihood at the mode
mode = x                                              # mode of the parameters
Sigma = - np.linalg.inv(aux.hessian(lik_f, mode)[0])  # inverse variance-covariance matrix

# MH algorithm
t = time.time()
para_sim, para_avg, para_5, para_95, logposterior, acceptancerate = aux.metropolis_hastings(lik_f, mode, Sigma, bounds,
                                                                                            Nsim, Nburn, c)
times['mh'] = time.time() - t
times['c'], times['rate'] = c, acceptancerate

# Compute recursive mean
rmean = np.zeros((para_avg.shape[0], para_sim.shape[1]))
for i in range(para_sim.shape[1]):
    rmean[:, i] = np.mean(para_sim[:, :i + 1], 1)

# Compute kernel density
kde = sm.nonparametric.KDEUnivariate(para_sim[0])
kde.fit()
kernel_support = np.zeros((para_avg.shape[0], kde.support.shape[0]))
kernel_density = np.zeros((para_avg.shape[0], kde.support.shape[0]))
for i in range(para_avg.shape[0]):
    kde = sm.nonparametric.KDEUnivariate(para_sim[i])
    kde.fit()
    kernel_support[i] = kde.support
    kernel_density[i] = kde.density

# Show times
t1 = times['mh']
print(f'Time for Metropolis-Hastings algorithm: {round(t1, 2)}s')
print(f'Acceptance rate: {round(acceptancerate, 2)}')
Time for Metropolis-Hastings algorithm: 61.5s
Acceptance rate: 0.24
In [31]:
# save data
if savedata: np.save('import_export/results/hank_1afull_est',
                     [x, para_avg, para_5, para_95, Nsim, acceptancerate, rmean, kernel_support, kernel_density,
                      para_sim])
if savedata: np.save('import_export/results/hank_1afull_speed', [times])

5. Two-asset HANK

5.1 Two-asset HANK (only shock processes)

Calibrate steady state and precompute Jacobian.

In [32]:
times = {}

# Steady state burn-in
t = time.time()
ss = hank_2a.hank_ss(noisy=False)
elapsed_ss0 = time.time() - t

# Steady state
t = time.time()
ss = hank_2a.hank_ss(noisy=False)
times['ss'] = time.time() - t

# Compute model jacobian
T = 300
exogenous = ['rstar', 'markup', 'G', 'beta', 'Z', 'rinv_shock', 'markup_w']
unknowns = ['r', 'w', 'Y']
targets = ['asset_mkt', 'fisher', 'wnkpc']
outputs = ['Y', 'pi', 'i', 'C', 'N', 'I', 'w']

t = time.time()
pricing = solved(block_list=[hank_2a.pricing], unknowns=['pi'], targets=['nkpc'])
arbitrage = solved(block_list=[hank_2a.arbitrage], unknowns=['p'], targets=['equity'])
production = solved(block_list=[hank_2a.labor, hank_2a.investment], unknowns=['Q', 'K'], targets=['inv', 'val'])
block_list = [hank_2a.household_inc, pricing, arbitrage, production, hank_2a.dividend, hank_2a.taylor, hank_2a.fiscal,
              hank_2a.finance, hank_2a.wage, hank_2a.union, hank_2a.mkt_clearing]
G = jac.get_G(block_list, exogenous, unknowns, targets, T=T, ss=ss, outputs=outputs)
elapsed_G0 = time.time() - t

t = time.time()
pricing = solved(block_list=[hank_2a.pricing], unknowns=['pi'], targets=['nkpc'])
arbitrage = solved(block_list=[hank_2a.arbitrage], unknowns=['p'], targets=['equity'])
production = solved(block_list=[hank_2a.labor, hank_2a.investment], unknowns=['Q', 'K'], targets=['inv', 'val'])
block_list = [hank_2a.household_inc, pricing, arbitrage, production, hank_2a.dividend, hank_2a.taylor, hank_2a.fiscal,
              hank_2a.finance, hank_2a.wage, hank_2a.union, hank_2a.mkt_clearing]
G = jac.get_G(block_list, exogenous, unknowns, targets, T=T, ss=ss, outputs=outputs)
times['G'] = time.time() - t

# Show times
t1, t2 = times['ss'], times['G']
print(f'\nTime to compute steady state (burn-in): {round(elapsed_ss0, 2)}s')
print(f'Time to compute model jacobian (burn-in): {round(elapsed_G0, 2)}s')

print(f'\nTime to compute steady state: {round(t1, 2)}s')
print(f'Time to compute model jacobian: {round(t2, 2)}s')
Time to compute steady state (burn-in): 17.75s
Time to compute model jacobian (burn-in): 5.24s

Time to compute steady state: 16.73s
Time to compute model jacobian: 4.75s

Nonlinear impulse response

In [33]:
rho_r = 0.6
drstar = -0.0125 * rho_r ** np.arange(T)  # to get pp change in rate, multiply by 4 since model is quarterly
td_nonlinear = nonlinear.td_solve(ss, block_list, unknowns, targets, rstar=ss['r'] + drstar, noisy=False)['C'] - ss['C']
td_linear = G['C']['rstar'] @ drstar

Transition to new steady state

In [34]:
# Compute new steady state
ssnew = hank_2a.hank_ss_Z(1.01, ss)

# Compute transition
td_newss = nonlinear.td_solve(ssnew, block_list, unknowns, targets, ss_initial=ss, Z=np.full(500, ssnew['Z']), noisy=False)

Bayesian estimation

In [35]:
# Define data
series = ['y', 'pi', 'i', 'c', 'n', 'I', 'w']
data = aux.get_normalized_data(ss, 'import_export/data/data_bayes.csv', series)

# Parameters of estimation
shock_series = [('Z', 1), ('rstar', 1), ('G', 1), ('beta', 1), ('rinv_shock', 1), ('markup', 1),
                ('markup_w', 1)]  # specifies shock to hit economy and whether it's an AR-1 or AR-2
x_guess = [0.4, 0.6] * len(shock_series)
priors = [('Invgamma', 0.4, 4), ('Beta', 0.5, 0.2)] * len(shock_series)
bounds = [(1e-2, 10), (1e-2, 0.98)] * len(shock_series)

# Estimation
t = time.time()
result, x_sd, times['nfev'] = aux.estimate(outputs, data, x_guess, shock_series, priors, T, G=G, bounds=bounds)
x = result.x
times_est0 = time.time() - t

t = time.time()
result, x_sd, times['nfev'] = aux.estimate(outputs, data, x_guess, shock_series, priors, T, G=G, bounds=bounds)
x = result.x
times['est'] = time.time() - t

# Show times
t1 = times['est']
print(f'Time for estimation (burn-in): {round(times_est0, 2)}s')
print(f'Time for estimation: {round(t1, 2)}s')
Time for estimation (burn-in): 14.98s
Time for estimation: 15.57s

Simulation at the posterior mode

In [36]:
T_simul = 100_000
t = time.time()
sim = aux.simulate(G, ['Z', 'rstar', 'G', 'beta', 'rinv_shock', 'markup', 'markup_w'], outputs,
                   {'Z': x[0], 'rstar': x[2], 'G': x[4], 'beta': x[6], 'rinv_shock': x[8], 'markup': x[10],
                    'markup_w': x[12]},
                   {'Z': x[1], 'rstar': x[3], 'G': x[5], 'beta': x[7], 'rinv_shock': x[9], 'markup': x[11],
                    'markup_w': x[13]}, T, T_simul)  # note, just simulate AR(1) for now
times['sim'] = time.time() - t

# Show times
t1 = times['sim']
print(f'Time for simulation: {round(1000 * t1, 3)}ms')
Time for simulation: 106.222ms

Get run time for likelihood evaluation

In [37]:
T_irf = T - 20
n_se,n_sh = len(outputs), len(shock_series)
meas_error = np.zeros(n_se)
sigmas, arcoefs, macoefs = aux.get_shocks_arma(x, shock_series)

# run them all once
As = aux.step1_est(G, arcoefs, macoefs, shock_series, outputs, T, T_irf, n_se, n_sh)
Sigma = ssj_est.all_covariances(As, sigmas)
llh = aux.step3_est(Sigma, data, sigma_o=meas_error)

# Compute run time
res = %timeit -o -q As = aux.step1_est(G, arcoefs, macoefs, shock_series, outputs, T, T_irf, n_se, n_sh)
times['lk_step1'] = np.mean(res.timings)
res = %timeit -o -q Sigma = ssj_est.all_covariances(As, sigmas)
times['lk_step2'] = np.mean(res.timings)
res = %timeit -o -q llh = aux.step3_est(Sigma, data, sigma_o=meas_error)
times['lk_step3'] = np.mean(res.timings)
times['lk'] = times['lk_step1'] + times['lk_step2'] + times['lk_step3']

Metropolis Hastings

In [ ]:
# Parameters
c = 0.4             # scale for the variance
x0 = x.copy()


# Define function that returns log-likelihood + log-prior
def lik_f(x):
    return aux.loglik_f(x, data, outputs, shock_series, priors, T, G)


# Get mode and variance covariance matrix of the log likelihood at the mode
mode = x                                              # mode of the parameters
Sigma = - np.linalg.inv(aux.hessian(lik_f, mode)[0])  # inverse variance-covariance matrix

# MH algorithm
t = time.time()
para_sim, para_avg, para_5, para_95, logposterior, acceptancerate = aux.metropolis_hastings(lik_f, mode, Sigma, bounds,
                                                                                            Nsim, Nburn, c)
times['mh'] = time.time() - t
times['c'], times['rate'] = c, acceptancerate

# Compute recursive mean
rmean = np.zeros((para_avg.shape[0], para_sim.shape[1]))
for i in range(para_sim.shape[1]):
    rmean[:, i] = np.mean(para_sim[:, :i + 1], 1)

# Compute kernel density
kde = sm.nonparametric.KDEUnivariate(para_sim[0])
kde.fit()
kernel_support = np.zeros((para_avg.shape[0], kde.support.shape[0]))
kernel_density = np.zeros((para_avg.shape[0], kde.support.shape[0]))
for i in range(para_avg.shape[0]):
    kde = sm.nonparametric.KDEUnivariate(para_sim[i])
    kde.fit()
    kernel_support[i] = kde.support
    kernel_density[i] = kde.density

# Show times
t1 = times['mh']
print(f'Time for Metropolis-Hastings algorithm: {round(t1, 2)}s')
print(f'Acceptance rate: {round(acceptancerate, 2)}')
In [38]:
# save data
if savedata: np.save('import_export/results/hank_2a_transition', [td_nonlinear, td_linear, td_newss, ss])
if savedata: np.save('import_export/results/hank_2a_est',
                     [x, para_avg, para_5, para_95, Nsim, acceptancerate, rmean, kernel_support, kernel_density,
                      para_sim])
if savedata: np.save('import_export/results/hank_2a_speed', [times])

5.2 Two-asset HANK model (shocks + parameters)

Bayesian estimation

In [39]:
times = {}

# Steady state
t = time.time()
ss = hank_2a.hank_ss(noisy=False)
times['ss'] = time.time() - t

# Define data
series = ['y', 'pi', 'i', 'c', 'n', 'I', 'w']
outputs = ['Y', 'pi', 'i', 'C', 'N', 'I', 'w']
data = aux.get_normalized_data(ss, 'import_export/data/data_bayes.csv', series)

# Options for estimation
params = ['phi', 'phi_y', 'kappap', 'kappaw', 'epsI']
shock_series = [('Z', 1), ('rstar', 1), ('G', 1), ('beta', 1), ('rinv_shock', 1), ('markup', 1),
                ('markup_w', 1)]  # specifies shock to hit economy and whether it's an AR-1 or AR-2
priors = [('Invgamma', 0.4, 4), ('Beta', 0.5, 0.2)] * len(shock_series) + [('Gamma', 1.5, 0.25), ('Gamma', 0.5, 0.25),
                                                                           ('Gamma', 0.1, 0.05), ('Gamma', 0.1, 0.05),
                                                                           ('Gamma', 4, 2)]
bounds = [(1e-2, 10), (1e-2, 0.99)] * len(shock_series) + [(1, 3), (1e-3, 4), (1e-4, 0.7), (1e-4, 0.7), (1e-3, 5)]
x_guess = [0.4, 0.6] * len(shock_series) + [ss[k] for k in params]
T = 300

# Compute jacobian of household block
J_ha = jac.get_G([hank_2a.household_inc], ['ra', 'rb', 'beta', 'w', 'N', 'tax'], [], [], T, ss)

# Store GE information
unknowns = ['r', 'w', 'Y']
targets = ['asset_mkt', 'fisher', 'wnkpc']
exogenous = ['rstar', 'markup', 'G', 'beta', 'Z', 'rinv_shock', 'markup_w']
pricing = solved(block_list=[hank_2a.pricing], unknowns=['pi'], targets=['nkpc'])
arbitrage = solved(block_list=[hank_2a.arbitrage], unknowns=['p'], targets=['equity'])
production = solved(block_list=[hank_2a.labor, hank_2a.investment], unknowns=['Q', 'K'], targets=['inv', 'val'])
block_list = [J_ha, pricing, arbitrage, production, hank_2a.dividend, hank_2a.taylor, hank_2a.fiscal, hank_2a.finance,
              hank_2a.wage, hank_2a.union, hank_2a.mkt_clearing]
jac_info = {'T': T, 'block_list': block_list, 'exogenous': exogenous, 'unknowns': unknowns, 'targets': targets, 'ss': ss}

# Estimate
t = time.time()
result, x_sd, times['nfev'] = aux.estimate(outputs, data, x_guess, shock_series, priors, T, params=params,
                                           jac_info=jac_info, bounds=bounds, sd=False)
x = result.x
times['est'] = time.time() - t

# Show times
t1, t2 = times['ss'], times['est']
print(f'\nTime to compute steady state: {round(t1, 2)}s')
print(f'Time for estimation: {round(t2, 2)}s')
Time to compute steady state: 17.56s
Time for estimation: 596.23s

Get run time for likelihood evaluation

In [40]:
n_se,n_sh = len(outputs), len(shock_series)
meas_error = np.zeros(n_se)
sigmas, arcoefs, macoefs = aux.get_shocks_arma(x, shock_series)

# run them all once
As = aux.step1_estfull(x,shock_series,outputs,T,params,jac_info)
Sigma = ssj_est.all_covariances(As, sigmas)
llh = aux.step3_est(Sigma, data, sigma_o=meas_error)

# Compute run time
res = %timeit -o -q As = aux.step1_estfull(x,shock_series,outputs,T,params,jac_info)
times['lk_step1'] = np.mean(res.timings)
res = %timeit -o -q Sigma = ssj_est.all_covariances(As, sigmas)
times['lk_step2'] = np.mean(res.timings)
res = %timeit -o -q llh = aux.step3_est(Sigma, data, sigma_o=meas_error)
times['lk_step3'] = np.mean(res.timings)
times['lk'] = times['lk_step1'] + times['lk_step2'] + times['lk_step3']

Metropolis Hastings

In [ ]:
# Parameters
c = 0.10            # scale for the variance
x0 = x.copy()

# Define function that returns log-likelihood + log-prior
def lik_f(x):
    n_params = len(params)
    x_params = x[-n_params:]

    # write new parameters into ss
    ss_here = jac_info['ss'].copy()
    ss_here.update({param: x_params[j] for j, param in enumerate(params)})

    # Compute model jacobian G
    G = jac.get_G(jac_info['block_list'], exogenous=jac_info['exogenous'], unknowns=jac_info['unknowns'],
                  targets=jac_info['targets'], T=jac_info['T'], ss=ss_here, outputs=outputs)

    # return likelihood
    return aux.loglik_f(x, data, outputs, shock_series, priors, T, G)


# Get mode and variance covariance matrix of the log likelihood at the mode
mode = x                                              # mode of the parameters
Sigma = - np.linalg.inv(aux.hessian(lik_f, mode)[0])  # inverse variance-covariance matrix

# MH algorithm
t = time.time()
para_sim, para_avg, para_5, para_95, logposterior, acceptancerate = aux.metropolis_hastings(lik_f, mode, Sigma, bounds,
                                                                                            Nsim, Nburn, c)
times['mh'] = time.time() - t
times['c'], times['rate'] = c, acceptancerate

# Compute recursive mean
rmean = np.zeros((para_avg.shape[0], para_sim.shape[1]))
for i in range(para_sim.shape[1]):
    rmean[:, i] = np.mean(para_sim[:, :i + 1], 1)

# Compute kernel density
kde = sm.nonparametric.KDEUnivariate(para_sim[0])
kde.fit()
kernel_support = np.zeros((para_avg.shape[0], kde.support.shape[0]))
kernel_density = np.zeros((para_avg.shape[0], kde.support.shape[0]))
for i in range(para_avg.shape[0]):
    kde = sm.nonparametric.KDEUnivariate(para_sim[i])
    kde.fit()
    kernel_support[i] = kde.support
    kernel_density[i] = kde.density

# Show times
t1 = times['mh']
print(f'Time for Metropolis-Hastings algorithm: {round(t1, 2)}s')
print(f'Acceptance rate: {round(acceptancerate, 2)}')
In [41]:
# save data
if savedata: np.save('import_export/results/hank_2afull_est',
                     [x, para_avg, para_5, para_95, Nsim, acceptancerate, rmean, kernel_support, kernel_density,
                      para_sim])
if savedata: np.save('import_export/results/hank_2afull_speed', [times])

6. Accuracy checks

6.1 Truncation error

In [42]:
# parameters
T_list = np.array([100, 125, 150, 175, 200, 300, 400, 500, 1000])
Ttrunc = 100
sigma = 0.01

# Compute stead states
ss_ks = ks.ks_ss()
ss_ha = hank_1a.hank_ss()
ss_ha2 = hank_2a.hank_ss(noisy=False)


def truncation_f(rho):
    # Compute model for different grid sizes
    irf_list_ks = np.zeros((Ttrunc, len(T_list)))
    irf_list_hank = np.zeros((Ttrunc, len(T_list)))
    irf_list_hank2 = np.zeros((Ttrunc, len(T_list)))
    for i in range(len(T_list)):
        # Krusell Smith
        G = jac.get_G(block_list=[ks.firm, ks.mkt_clearing, ks.household], exogenous=['Z'], unknowns=['K'],
                      targets=['asset_mkt'], T=T_list[i], ss=ss_ks)
        dZ = ss_ks['Z'] * sigma * rho ** np.arange(T_list[i])
        irf_list_ks[:, i] = (G['Y']['Z'] @ dZ / ss_ks['Y'])[:Ttrunc]

        # one asset hank
        G = jac.get_G([hank_1a.household_trans, hank_1a.firm, hank_1a.monetary, hank_1a.fiscal, hank_1a.nkpc,
                       hank_1a.mkt_clearing],
                      exogenous=['Z'], unknowns=['pi', 'w', 'Y'], targets=['nkpc_res', 'asset_mkt', 'labor_mkt'],
                      T=T_list[i], ss=ss_ha)
        dZ = ss_ha['Z'] * sigma * rho ** np.arange(T_list[i])
        irf_list_hank[:, i] = (G['Y']['Z'] @ dZ / ss_ha['Y'])[:Ttrunc]

        # two asset hank
        pricing = solved(block_list=[hank_2a.pricing], unknowns=['pi'], targets=['nkpc'])
        arbitrage = solved(block_list=[hank_2a.arbitrage], unknowns=['p'], targets=['equity'])
        production = solved(block_list=[hank_2a.labor, hank_2a.investment], unknowns=['Q', 'K'], targets=['inv', 'val'])
        block_list = [hank_2a.household_inc, pricing, arbitrage, production, hank_2a.dividend, hank_2a.taylor,
                      hank_2a.fiscal, hank_2a.finance, hank_2a.wage, hank_2a.union, hank_2a.mkt_clearing]
        G = jac.get_G(block_list, exogenous=['rstar', 'markup', 'G', 'beta', 'Z', 'rinv_shock', 'markup_w'],
                      unknowns=['r', 'w', 'Y'], targets=['asset_mkt', 'fisher', 'wnkpc'], T=T_list[i], ss=ss_ha2,
                      outputs=['Y'])
        dZ = ss_ha2['Z'] * sigma * rho ** np.arange(T_list[i])
        irf_list_hank2[:, i] = (G['Y']['Z'] @ dZ / ss_ha2['Y'])[:Ttrunc]

    # Compute mean squared error
    diff_ks = irf_list_ks - np.tile(irf_list_ks[:, -1], (len(T_list), 1)).T
    diff_hank = irf_list_hank - np.tile(irf_list_hank[:, -1], (len(T_list), 1)).T
    diff_hank2 = irf_list_hank2 - np.tile(irf_list_hank2[:, -1], (len(T_list), 1)).T
    rmse_ks = 1 / sigma * np.sqrt(np.sum(diff_ks ** 2, axis=0) / diff_ks.shape[0])  # Report in percent
    rmse_hank = 1 / sigma * np.sqrt(np.sum(diff_hank ** 2, axis=0) / diff_hank.shape[0])  # Report in percent
    rmse_hank2 = 1 / sigma * np.sqrt(np.sum(diff_hank2 ** 2, axis=0) / diff_hank2.shape[0])  # Report in percent

    return rmse_ks, rmse_hank, rmse_hank2


# Compute error for rho = 0.9 and rho = 0.90
rmse_ks_rho90, rmse_hank_rho90, rmse_hank2_rho90 = truncation_f(0.9)
rmse_ks_rho99, rmse_hank_rho99, rmse_hank2_rho99 = truncation_f(0.99)

# save data
if savedata: np.save('import_export/results/truncation',
                     [T_list, rmse_ks_rho90, rmse_hank_rho90, rmse_hank2_rho90, rmse_ks_rho99, rmse_hank_rho99,
                      rmse_hank2_rho99])

6.2 Equivalence between fake news algorithm and direct method

In [43]:
# Steady state
ss_ks = ks.ks_ss()
ss_ha = hank_1a.hank_ss()
ss_ha2 = hank_2a.hank_ss(noisy=False)

# Parameters for irf
cols = np.array([5, 25, 50, 100])
T = 150

# Compute jacobian with fake news algorithm
J_ks = jac.get_G([ks.household], ['r', 'w'], [], [], T, ss_ks)
J_ha = jac.get_G([hank_1a.household_trans], ['r', 'w', 'Div', 'Tax'], [], [], T, ss_ha)
J_ha2 = jac.get_G([hank_2a.household_inc], ['ra', 'rb', 'beta', 'w', 'N', 'tax'], [], [], T, ss_ha2)

# Compute jacobian directly
J_ks_direct = aux_jac.get_J_direct(ks.household, ['r', 'w'], ['C'], ss_ks, T, cols)
J_ha_direct = aux_jac.get_J_direct(hank_1a.household_trans, ['r', 'w', 'Div', 'Tax'], ['C'], ss_ha, T, cols)
J_ha2_direct = aux_jac.get_J_direct(hank_2a.household_inc, ['ra', 'rb', 'beta', 'w', 'N', 'tax'], ['C'], ss_ha2, T, cols)

# save data
if savedata: np.save('import_export/results/fake_news',
                     [cols, J_ks, J_ks_direct, J_ha, J_ha_direct, J_ha2, J_ha2_direct])

7 Computing times

7.1. HA Jacobian

In [44]:
times = {'direct': {0: {}, 1: {}, 2: {}}, 'fake': {0: {}, 1: {}, 2: {}, 3: {}, 4: {}}}

Krusell-Smith

$n_x=2, n_y=2, n_g=3500.$

In [45]:
# Steady state
ss_ks = ks.ks_ss()
T = 300

## Fake-news algorithm

# Call all the steps once
shock_dict = {'r': {'r': 1}, 'w': {'w': 1}}
ssinput_dict, ssy_list, outcome_list, V_name, a_pol_i, a_pol_pi, a_space = aux_jac.step_fake_0(ks.backward_iterate, ss_ks)
curlyYs, curlyDs = aux_jac.step_fake_1(ks.backward_iterate, {'r': {'r': 1}, 'w': {'w': 1}}, ss_ks, ssinput_dict, ssy_list, outcome_list, V_name, a_pol_i, a_space, T)
curlyPs = aux_jac.step_fake_2(ss_ks, outcome_list, ssy_list, a_pol_i, a_pol_pi, T)
F = aux_jac.step_fake_3(shock_dict, outcome_list, curlyYs, curlyDs, curlyPs)
J = aux_jac.step_fake_4(F, shock_dict, outcome_list)

# Step 1
res = %timeit -o -q curlyYs, curlyDs = aux_jac.step_fake_1(ks.backward_iterate, {'r': {'r': 1}, 'w': {'w': 1}}, ss_ks, ssinput_dict, ssy_list, outcome_list, V_name, a_pol_i, a_space, T)
times['fake'][1]['ks'] = np.mean(res.timings)

# Step 2
res = %timeit -o -q curlyPs = aux_jac.step_fake_2(ss_ks, outcome_list, ssy_list, a_pol_i, a_pol_pi, T)
times['fake'][2]['ks'] = np.mean(res.timings)

# Step 3
res = %timeit -o -q F = aux_jac.step_fake_3(shock_dict, outcome_list, curlyYs, curlyDs, curlyPs)
times['fake'][3]['ks'] = np.mean(res.timings)

# Step 4
res = %timeit -o -q J = aux_jac.step_fake_4(F, shock_dict, outcome_list)
times['fake'][4]['ks'] = np.mean(res.timings)

# Compute total
avg = 0
for i in range(1,5): avg += np.mean(times['fake'][i]['ks'])
times['fake'][0]['ks'] = avg 

## Direct method

# Call all the steps once
Va_path, a_path, c_path, D_path = aux_jac.ks_backward(ss_ks, ks.backward_iterate, w=np.full(T, ss_ks['w']) + 1E-4*(np.arange(T) == 3)) 
col1 = aux_jac.ks_forward(ss_ks, T, Va_path, a_path, c_path, D_path)

# Step 1
res = %timeit -o -q Va_path, a_path, c_path, D_path = aux_jac.ks_backward(ss_ks, ks.backward_iterate, w=np.full(T, ss_ks['w']) + 1E-4*(np.arange(T) == 3)) 
times['direct'][1]['ks'] = np.mean(res.timings) * 600 
# multiply by 600 because there are 600 jacobians to compute (2 inputs x 300 periods)

# Step 2
res = %timeit -o -q col1 = aux_jac.ks_forward(ss_ks, T, Va_path, a_path, c_path, D_path)
times['direct'][2]['ks'] = np.mean(res.timings) * 600
# multiply by 600 because there are 600 jacobians to compute (2 inputs x 300 periods)

# Compute total
avg = 0
for i in range(1,3): avg += np.mean(times['direct'][i]['ks'])
times['direct'][0]['ks'] = avg

High-dimensional Krusell-Smith

$n_x=2, n_y=2, n_g=250000.$

In [46]:
# Steady state
ss_ks = ks.ks_ss(nS=50, nA=5000)
T = 300

## Fake-news algorithm

# Call all the steps once
shock_dict = {'r': {'r': 1}, 'w': {'w': 1}}
ssinput_dict, ssy_list, outcome_list, V_name, a_pol_i, a_pol_pi, a_space = aux_jac.step_fake_0(ks.backward_iterate, ss_ks)
curlyYs, curlyDs = aux_jac.step_fake_1(ks.backward_iterate, {'r': {'r': 1}, 'w': {'w': 1}}, ss_ks, ssinput_dict, ssy_list,
                                       outcome_list, V_name, a_pol_i, a_space, T)
curlyPs = aux_jac.step_fake_2(ss_ks, outcome_list, ssy_list, a_pol_i, a_pol_pi, T)
F = aux_jac.step_fake_3(shock_dict, outcome_list, curlyYs, curlyDs, curlyPs)
J = aux_jac.step_fake_4(F, shock_dict, outcome_list)

# Step 1
res = %timeit -o -q curlyYs, curlyDs = aux_jac.step_fake_1(ks.backward_iterate, {'r': {'r': 1}, 'w': {'w': 1}}, ss_ks, ssinput_dict, ssy_list, outcome_list, V_name, a_pol_i, a_space, T)
times['fake'][1]['hd_ks'] = np.mean(res.timings)

# Step 2
res = %timeit -o -q curlyPs = aux_jac.step_fake_2(ss_ks, outcome_list, ssy_list, a_pol_i, a_pol_pi, T)
times['fake'][2]['hd_ks'] = np.mean(res.timings)

# Step 3
res = %timeit -o -q F = aux_jac.step_fake_3(shock_dict, outcome_list, curlyYs, curlyDs, curlyPs)
times['fake'][3]['hd_ks'] = np.mean(res.timings)

# Step 4
res = %timeit -o -q J = aux_jac.step_fake_4(F, shock_dict, outcome_list)
times['fake'][4]['hd_ks'] = np.mean(res.timings)

# Compute total
avg = 0
for i in range(1,5): avg += np.mean(times['fake'][i]['hd_ks'])
times['fake'][0]['hd_ks'] = avg 

## Direct method

# Call all the steps once
Va_path, a_path, c_path, D_path = aux_jac.ks_backward(ss_ks, ks.backward_iterate, w=np.full(T, ss_ks['w']) + 1E-4*(np.arange(T) == 3)) 
col1 = aux_jac.ks_forward(ss_ks, T, Va_path, a_path, c_path, D_path)

# Step 1
res = %timeit -o -q Va_path, a_path, c_path, D_path = aux_jac.ks_backward(ss_ks, ks.backward_iterate, w=np.full(T, ss_ks['w']) + 1E-4*(np.arange(T) == 3)) 
times['direct'][1]['hd_ks'] = np.mean(res.timings) * 600
# multiply by 600 because there are 600 jacobians to compute (2 inputs x 300 periods)

# Step 2
res = %timeit -o -q col1 = aux_jac.ks_forward(ss_ks, T, Va_path, a_path, c_path, D_path)
times['direct'][2]['hd_ks'] = np.mean(res.timings) * 600
# multiply by 600 because there are 600 jacobians to compute (2 inputs x 300 periods)

# Compute total
avg = 0
for i in range(1,3): avg += np.mean(times['direct'][i]['hd_ks'])
times['direct'][0]['hd_ks'] = avg

One-asset hank

$n_x=4, n_y=4, n_g=3500.$

In [47]:
# Steady state
ss_ha = hank_1a.hank_ss()
T = 300

## Fake-news algorithm

# Call all the steps once
T_div = ss_ha['div_rule'] / np.sum(ss_ha['pi_e'] * ss_ha['div_rule'])
T_tax = -ss_ha['tax_rule'] / np.sum(ss_ha['pi_e'] * ss_ha['tax_rule'])
shock_dict = {'r': {'r': 1}, 'w': {'w': 1}, 'Div': {'T': T_div}, 'Tax': {'T': T_tax}}
ssinput_dict, ssy_list, outcome_list, V_name, a_pol_i, a_pol_pi, a_space = aux_jac.step_fake_0(hank_1a.backward_iterate, ss_ha)
curlyYs, curlyDs = aux_jac.step_fake_1(hank_1a.backward_iterate, shock_dict, ss_ha, ssinput_dict, ssy_list, outcome_list, V_name, a_pol_i, a_space, T)
curlyPs = aux_jac.step_fake_2(ss_ha, outcome_list, ssy_list, a_pol_i, a_pol_pi, T)
F = aux_jac.step_fake_3(shock_dict, outcome_list, curlyYs, curlyDs, curlyPs)
J = aux_jac.step_fake_4(F, shock_dict, outcome_list)

# Step 1
res = %timeit -o -q curlyYs, curlyDs = aux_jac.step_fake_1(hank_1a.backward_iterate, {'r': {'r': 1}, 'w': {'w': 1}}, ss_ha, ssinput_dict, ssy_list, outcome_list, V_name, a_pol_i, a_space, T)
times['fake'][1]['ha'] = np.mean(res.timings)

# Step 2
res = %timeit -o -q curlyPs = aux_jac.step_fake_2(ss_ha, outcome_list, ssy_list, a_pol_i, a_pol_pi, T)
times['fake'][2]['ha'] = np.mean(res.timings)

# Step 3
res = %timeit -o -q F = aux_jac.step_fake_3(shock_dict, outcome_list, curlyYs, curlyDs, curlyPs)
times['fake'][3]['ha'] = np.mean(res.timings)

# Step 4
res = %timeit -o -q J = aux_jac.step_fake_4(F, shock_dict, outcome_list)
times['fake'][4]['ha'] = np.mean(res.timings)

# Compute total
avg = 0
for i in range(1,5): avg += np.mean(times['fake'][i]['ha'])
times['fake'][0]['ha'] = avg 

## Direct method

# Call all the steps once
Va_path, a_path, c_path, n_path, ns_path, D_path = aux_jac.ha_backward(ss_ha, hank_1a.backward_iterate, w=np.full(T, ss_ha['w']) + 1E-4*(np.arange(T)==3))
col1 = aux_jac.ha_forward(ss_ha, T, Va_path, a_path, c_path, n_path, ns_path, D_path)

# Step 1
res = %timeit -o -q Va_path, a_path, c_path, n_path, ns_path, D_path = aux_jac.ha_backward(ss_ha, hank_1a.backward_iterate, w=np.full(T, ss_ha['w']) + 1E-4*(np.arange(T)==3))
times['direct'][1]['ha'] = np.mean(res.timings) * 1200
# multiply by 1200 because there are 1200 jacobians to compute (4 inputs x 300 periods)

# Step 2
res = %timeit -o -q col1 = aux_jac.ha_forward(ss_ha, T, Va_path, a_path, c_path, n_path, ns_path, D_path)
times['direct'][2]['ha'] = np.mean(res.timings) * 1200
# multiply by 1200 because there are 1200 jacobians to compute (4 inputs x 300 periods)

# Compute total
avg = 0
for i in range(1,3): avg += np.mean(times['direct'][i]['ha'])
times['direct'][0]['ha'] = avg

Two-asset hank

$n_x=5, n_y=4, n_g=10500$

In [48]:
# Compute stead states
ss_ha2 = hank_2a.hank_ss(noisy=False)
T = 300

## Fake-news algorithm

# Call all the steps once
shock_dict = {'ra': {'ra': 1},
              'rb': {'rb': 1},
              'tax': {'z_grid': -ss_ha2['w']*ss_ha2['e_grid']},
              'w': {'z_grid': (1-ss_ha2['tax'])*ss_ha2['e_grid']},
              'N': {'z_grid': (1-ss_ha2['tax'])*ss_ha2['w']*ss_ha2['e_grid']}}
ssinput_dict, ssy_list, outcome_list, V_list, a_i, a_pi, b_i, b_pi, a_space, b_space = aux_jac.step_fake_0_2d(hank_2a.backward_iterate, ss_ha2)
curlyYs, curlyDs = aux_jac.step_fake_1_2d(hank_2a.backward_iterate, shock_dict, ss_ha2, ssinput_dict, ssy_list, outcome_list, V_list, a_i, a_pi, b_i, b_pi, a_space, b_space, T)
curlyPs = aux_jac.step_fake_2_2d(ss_ha2, outcome_list, ssy_list, a_i, b_i, a_pi, b_pi, T)
F = aux_jac.step_fake_3_2d(shock_dict, outcome_list, curlyYs, curlyDs, curlyPs)
J = aux_jac.step_fake_4_2d(F, shock_dict, outcome_list)

# Step 1
res = %timeit -o -q curlyYs, curlyDs = aux_jac.step_fake_1_2d(hank_2a.backward_iterate, shock_dict, ss_ha2, ssinput_dict, ssy_list, outcome_list, V_list, a_i, a_pi, b_i, b_pi, a_space, b_space, T)
times['fake'][1]['ha2'] = np.mean(res.timings)

# Step 2
res = %timeit -o -q curlyPs = aux_jac.step_fake_2_2d(ss_ha2, outcome_list, ssy_list, a_i, b_i, a_pi, b_pi, T)
times['fake'][2]['ha2'] = np.mean(res.timings)

# Step 3
res = %timeit -o -q F = aux_jac.step_fake_3_2d(shock_dict, outcome_list, curlyYs, curlyDs, curlyPs)
times['fake'][3]['ha2'] = np.mean(res.timings)

# Step 4
res = %timeit -o -q J = aux_jac.step_fake_4_2d(F, shock_dict, outcome_list)
times['fake'][4]['ha2'] = np.mean(res.timings)

# Compute total
avg = 0
for i in range(1,5): avg += np.mean(times['fake'][i]['ha2'])
times['fake'][0]['ha2'] = avg 

## Direct method

# Call all the steps once
Va_t, Vb_t, a_t, b_t, c_t, m_t, D_treated, D_control = aux_jac.ha2_backward(ss_ha2, hank_2a.backward_iterate, ra=np.full(T, ss_ha2['ra']) + 1E-4*(np.arange(T)==3))
col3 = aux_jac.ha2_forward(ss_ha2, T, Va_t, Vb_t, a_t, b_t, c_t, m_t, D_treated, D_control)

# Step 1
res = %timeit -o -q Va_t, Vb_t, a_t, b_t, c_t, m_t, D_treated, D_control = aux_jac.ha2_backward(ss_ha2, hank_2a.backward_iterate, ra=np.full(T, ss_ha2['ra']) + 1E-4*(np.arange(T)==3))
times['direct'][1]['ha2'] = np.mean(res.timings) * 1500
# multiply by 1500 because there are 1500 jacobians to compute (5 inputs x 300 periods)

# Step 2
res = %timeit -o -q col3 = aux_jac.ha2_forward(ss_ha2, T, Va_t, Vb_t, a_t, b_t, c_t, m_t, D_treated, D_control)
times['direct'][2]['ha2'] = np.mean(res.timings) * 1500/2
# dividing by two because this method duplicates effort  

# Compute total
avg = 0
for i in range(1,3): avg += np.mean(times['direct'][i]['ha2'])
times['direct'][0]['ha2'] = avg
In [49]:
if savedata: np.save('import_export/results/fake_news_speed', [times])

7.2 General equilibrium Jacobian

In [50]:
T = 300
times = {'DAG': {0: {}, 1: {}, 2: {}, 3: {}}, 'no_DAG': {0: {}, 1: {}, 2: {}, 3: {}}}

Krusell-Smith

In [51]:
# Steady state
ss = ks.ks_ss()

# Compute jacobian of household block
J_ha = jac.get_G([ks.household], ['r', 'w'], [], [], T, ss)

## Using DAG

# define model
exogenous = ['Z']
unknowns = ['K']
targets = ['asset_mkt']
outputs = ['K', 'w', 'r', 'C', 'Y']
block_list=[J_ha, ks.firm, ks.mkt_clearing]

# run once
curlyJs, required = aux_jac.step0_preliminaries(block_list, exogenous, unknowns, ss)
J_curlyH = aux_jac.step1_forward_H(curlyJs, exogenous, unknowns, targets, required)
G_U = aux_jac.step2_solve_GU(J_curlyH, exogenous, unknowns, targets, T)
G = aux_jac.step3_forward_G(G_U, curlyJs, exogenous, targets, outputs, required, unknowns)

# Step 1
res = %timeit -o -q J_curlyH = aux_jac.step1_forward_H(curlyJs, exogenous, unknowns, targets, required)
times['DAG'][1]['ks'] = np.mean(res.timings)

# Step 2
res = %timeit -o -q G_U = aux_jac.step2_solve_GU(J_curlyH, exogenous, unknowns, targets, T)
times['DAG'][2]['ks'] = np.mean(res.timings)

# Step 3
res = %timeit -o -q G = aux_jac.step3_forward_G(G_U, curlyJs, exogenous, targets, outputs, required, unknowns)
times['DAG'][3]['ks'] = np.mean(res.timings)

# Compute total
avg = 0
for i in range(1,4): avg += np.mean(times['DAG'][i]['ks'])
times['DAG'][0]['ks'] = avg

## No DAG

# define model
unknowns = ['Y', 'r', 'w', 'K']
targets = ['res_Y', 'res_r', 'res_w', 'asset_mkt']
block_list = [J_ha, ks.ks_nodag]

# run once
curlyJs, required = aux_jac.step0_preliminaries(block_list, exogenous, unknowns, ss)
J_curlyH = aux_jac.step1_forward_H(curlyJs, exogenous, unknowns, targets, required)
G_U = aux_jac.step2_solve_GU(J_curlyH, exogenous, unknowns, targets, T)
G = aux_jac.step3_forward_G(G_U, curlyJs, exogenous, targets, outputs, required, unknowns)

# Step 1
res = %timeit -o -q J_curlyH = aux_jac.step1_forward_H(curlyJs, exogenous, unknowns, targets, required)
times['no_DAG'][1]['ks'] = np.mean(res.timings)

# Step 2
res = %timeit -o -q G_U = aux_jac.step2_solve_GU(J_curlyH, exogenous, unknowns, targets, T)
times['no_DAG'][2]['ks'] = np.mean(res.timings)

# Step 3
res = %timeit -o -q G = aux_jac.step3_forward_G(G_U, curlyJs, exogenous, targets, outputs, required, unknowns)
times['no_DAG'][3]['ks'] = np.mean(res.timings)

# Compute total
avg = 0
for i in range(1,4): avg += np.mean(times['no_DAG'][i]['ks'])
times['no_DAG'][0]['ks'] = avg

One-asset hank

In [52]:
# Steady state
ss = hank_1a.hank_ss()

# Compute jacobian of household block
J_ha = jac.get_G([hank_1a.household_trans], ['r', 'w', 'Div', 'Tax'], [], [], T, ss, save=True)

## Using DAG

# define model
exogenous = ['G', 'markup', 'rstar']
unknowns = ['pi', 'w', 'Y']
targets = ['nkpc_res', 'asset_mkt', 'labor_mkt']
outputs = ['pi', 'Tax', 'Y', 'Div', 'r', 'i', 'N', 'w', 'A', 'NS']
block_list=[J_ha, hank_1a.firm, hank_1a.monetary, hank_1a.fiscal, hank_1a.nkpc, hank_1a.mkt_clearing]

# run once
curlyJs, required = aux_jac.step0_preliminaries(block_list, exogenous, unknowns, ss)
J_curlyH = aux_jac.step1_forward_H(curlyJs, exogenous, unknowns, targets, required)
G_U = aux_jac.step2_solve_GU(J_curlyH, exogenous, unknowns, targets, T)
G = aux_jac.step3_forward_G(G_U, curlyJs, exogenous, targets, outputs, required, unknowns)

# Step 1
res = %timeit -o -q J_curlyH = aux_jac.step1_forward_H(curlyJs, exogenous, unknowns, targets, required)
times['DAG'][1]['ha'] = np.mean(res.timings)

# Step 2
res = %timeit -o -q G_U = aux_jac.step2_solve_GU(J_curlyH, exogenous, unknowns, targets, T)
times['DAG'][2]['ha'] = np.mean(res.timings)

# Step 3
res = %timeit -o -q G = aux_jac.step3_forward_G(G_U, curlyJs, exogenous, targets, outputs, required, unknowns)
times['DAG'][3]['ha'] = np.mean(res.timings)

# Compute total
avg = 0
for i in range(1,4): avg += np.mean(times['DAG'][i]['ha'])
times['DAG'][0]['ha'] = avg

## No DAG

# define model
unknowns = ['Y', 'L', 'w', 'pi', 'Div', 'Tax', 'r']
targets = ['res_L', 'res_Div', 'res_r', 'res_Tax', 'nkpc_res', 'asset_mkt', 'labor_mkt']
block_list=[J_ha, hank_1a.ha1_nodag]

# run once
curlyJs, required = aux_jac.step0_preliminaries(block_list, exogenous, unknowns, ss)
J_curlyH = aux_jac.step1_forward_H(curlyJs, exogenous, unknowns, targets, required)
G_U = aux_jac.step2_solve_GU(J_curlyH, exogenous, unknowns, targets, T)
G = aux_jac.step3_forward_G(G_U, curlyJs, exogenous, targets, outputs, required, unknowns)

# Step 1
res = %timeit -o -q J_curlyH = aux_jac.step1_forward_H(curlyJs, exogenous, unknowns, targets, required)
times['no_DAG'][1]['ha'] = np.mean(res.timings)

# Step 2
res = %timeit -o -q G_U = aux_jac.step2_solve_GU(J_curlyH, exogenous, unknowns, targets, T)
times['no_DAG'][2]['ha'] = np.mean(res.timings)

# Step 3
res = %timeit -o -q G = aux_jac.step3_forward_G(G_U, curlyJs, exogenous, targets, outputs, required, unknowns)
times['no_DAG'][3]['ha'] = np.mean(res.timings)

# Compute total
avg = 0
for i in range(1,4): avg += np.mean(times['no_DAG'][i]['ha'])
times['no_DAG'][0]['ha'] = avg

Two-asset hank

In [53]:
# Steady state
ss = hank_2a.hank_ss(noisy=False)

# Compute jacobian of household block
J_ha = jac.get_G([hank_2a.household_inc], ['ra', 'rb', 'beta', 'w', 'N', 'tax'], [], [], T, ss)

## Using DAG

# define model
exogenous = ['rstar', 'markup', 'G', 'beta', 'Z', 'rinv_shock', 'markup_w']
unknowns = ['r', 'w', 'Y']
targets = ['asset_mkt', 'fisher', 'wnkpc']
G_pricing = jac.get_G([hank_2a.pricing], ['mc', 'r', 'Y'], ['pi'], ['nkpc'],  T, ss)
G_arbitrage = jac.get_G([hank_2a.arbitrage], ['div', 'r'], ['p'],['equity'], T, ss)
G_production = jac.get_G([hank_2a.labor, hank_2a.investment], ['Y', 'w', 'r', 'Z'], ['Q', 'K'], ['inv', 'val'], T, ss)
block_list = [J_ha, G_pricing, G_arbitrage, G_production, hank_2a.dividend, hank_2a.taylor, hank_2a.fiscal, hank_2a.finance, hank_2a.wage, hank_2a.union, hank_2a.mkt_clearing]
outputs = None

# run once
curlyJs, required = aux_jac.step0_preliminaries(block_list, exogenous, unknowns, ss)
J_curlyH = aux_jac.step1_forward_H(curlyJs, exogenous, unknowns, targets, required)
G_U = aux_jac.step2_solve_GU(J_curlyH, exogenous, unknowns, targets, T)
G = aux_jac.step3_forward_G(G_U, curlyJs, exogenous, targets, outputs, required, unknowns)

# Step 1
res = %timeit -o -q J_curlyH = aux_jac.step1_forward_H(curlyJs, exogenous, unknowns, targets, required)
times['DAG'][1]['ha2'] = np.mean(res.timings)

# Step 2
res = %timeit -o -q G_U = aux_jac.step2_solve_GU(J_curlyH, exogenous, unknowns, targets, T)
times['DAG'][2]['ha2'] = np.mean(res.timings)

# Step 3
res = %timeit -o -q G = aux_jac.step3_forward_G(G_U, curlyJs, exogenous, targets, outputs, required, unknowns)
times['DAG'][3]['ha2'] = np.mean(res.timings)

# Compute total
avg = 0
for i in range(1,4): avg += np.mean(times['DAG'][i]['ha2'])
times['DAG'][0]['ha2'] = avg

## No DAG

# define model
unknowns = ['piw', 'psiw', 'rb', 'ra', 'tax', 'i', 'psip', 'I', 'Q', 'w', 'N', 'K', 'div', 'p', 'pi', 'mc', 'r', 'Y']
targets = ['nkpc', 'equity', 'res_N', 'res_mc', 'inv', 'val', 'res_psip', 'res_I', 'res_div', 'res_i', 'res_tax', 'res_rb',
           'res_ra', 'fisher', 'res_piw', 'res_psiw', 'wnkpc', 'asset_mkt']
block_list = [J_ha, hank_2a.ha2_nodag]

# run once
curlyJs, required = aux_jac.step0_preliminaries(block_list, exogenous, unknowns, ss)
J_curlyH = aux_jac.step1_forward_H(curlyJs, exogenous, unknowns, targets, required)
G_U = aux_jac.step2_solve_GU(J_curlyH, exogenous, unknowns, targets, T)
G = aux_jac.step3_forward_G(G_U, curlyJs, exogenous, targets, outputs, required, unknowns)

# Step 1
res = %timeit -o -q J_curlyH = aux_jac.step1_forward_H(curlyJs, exogenous, unknowns, targets, required)
times['no_DAG'][1]['ha2'] = np.mean(res.timings)

# Step 2
res = %timeit -o -q G_U = aux_jac.step2_solve_GU(J_curlyH, exogenous, unknowns, targets, T)
times['no_DAG'][2]['ha2'] = np.mean(res.timings)

# Step 3
res = %timeit -o -q G = aux_jac.step3_forward_G(G_U, curlyJs, exogenous, targets, outputs, required, unknowns)
times['no_DAG'][3]['ha2'] = np.mean(res.timings)

# Compute total
avg = 0
for i in range(1,4): avg += np.mean(times['no_DAG'][i]['ha2'])
times['no_DAG'][0]['ha2'] = avg
In [54]:
if savedata: np.save('import_export/results/G_speed',[times])

7.3 Linearized impulse responses

In [55]:
T = 300
times = {'DAG': {0: {}, 1: {}, 2:{}, 3:{}, 4:{}}, 'no_DAG': {0: {}, 1: {}, 2:{}, 3:{}, 4:{}}}

Krusell-Smith

In [56]:
# Steady state
ss = ks.ks_ss()

# Compute jacobian of household block
J_ha = jac.get_G([ks.household], ['r', 'w'], [], [], T, ss)

## Using DAG

# define model
exogenous = ['Z']
unknowns = ['K']
targets = ['asset_mkt']
outputs = ['K', 'w', 'r', 'C', 'Y']
block_list=[J_ha, ks.firm, ks.mkt_clearing]
dZ = {'Z': 0.9**np.arange(T)}

# run once
curlyJs, required = aux_jac.step0_irf_preliminaries(block_list, dZ, unknowns, ss)
H_U_unpacked = aux_jac.step1_irf_forward_HU(curlyJs, unknowns, targets, required)
J_curlyZ_dZ = aux_jac.step2_irf_forward_Z(curlyJs, dZ, targets, outputs, required)
dU = aux_jac.step3_irf_solve_dU(H_U_unpacked, J_curlyZ_dZ, unknowns, targets, T)
dX = aux_jac.step4_irf_forward_U_combine(curlyJs, dU, J_curlyZ_dZ, outputs, required, T)

# Step 1
res = %timeit -o -q H_U_unpacked = aux_jac.step1_irf_forward_HU(curlyJs, unknowns, targets, required)
times['DAG'][1]['ks'] = np.mean(res.timings)

# Step 2
res = %timeit -o -q J_curlyZ_dZ = aux_jac.step2_irf_forward_Z(curlyJs, dZ, targets, outputs, required)
times['DAG'][2]['ks'] = np.mean(res.timings)

# Step 3
res = %timeit -o -q dU = aux_jac.step3_irf_solve_dU(H_U_unpacked, J_curlyZ_dZ, unknowns, targets, T)
times['DAG'][3]['ks'] = np.mean(res.timings)

# Step 4
res = %timeit -o -q dX = aux_jac.step4_irf_forward_U_combine(curlyJs, dU, J_curlyZ_dZ, outputs, required, T)
times['DAG'][4]['ks'] = np.mean(res.timings)

# Compute total
avg = 0
for i in range(1,5): avg += np.mean(times['DAG'][i]['ks'])
times['DAG'][0]['ks'] = avg

## No DAG

# define model
unknowns = ['Y', 'r', 'w', 'K']
targets = ['res_Y', 'res_r', 'res_w', 'asset_mkt']
block_list = [J_ha, ks.ks_nodag]

# run once
curlyJs, required = aux_jac.step0_irf_preliminaries(block_list, dZ, unknowns, ss)
H_U_unpacked = aux_jac.step1_irf_forward_HU(curlyJs, unknowns, targets, required)
J_curlyZ_dZ = aux_jac.step2_irf_forward_Z(curlyJs, dZ, targets, outputs, required)
dU = aux_jac.step3_irf_solve_dU(H_U_unpacked, J_curlyZ_dZ, unknowns, targets, T)
dX = aux_jac.step4_irf_forward_U_combine(curlyJs, dU, J_curlyZ_dZ, outputs, required, T)

# Step 1
res = %timeit -o -q H_U_unpacked = aux_jac.step1_irf_forward_HU(curlyJs, unknowns, targets, required)
times['no_DAG'][1]['ks'] = np.mean(res.timings)

# Step 2
res = %timeit -o -q J_curlyZ_dZ = aux_jac.step2_irf_forward_Z(curlyJs, dZ, targets, outputs, required)
times['no_DAG'][2]['ks'] = np.mean(res.timings)

# Step 3
res = %timeit -o -q dU = aux_jac.step3_irf_solve_dU(H_U_unpacked, J_curlyZ_dZ, unknowns, targets, T)
times['no_DAG'][3]['ks'] = np.mean(res.timings)

# Step 4
res = %timeit -o -q dX = aux_jac.step4_irf_forward_U_combine(curlyJs, dU, J_curlyZ_dZ, outputs, required, T)
times['no_DAG'][4]['ks'] = np.mean(res.timings)

# Compute total
avg = 0
for i in range(1,5): avg += np.mean(times['no_DAG'][i]['ks'])
times['no_DAG'][0]['ks'] = avg

One-asset hank

In [57]:
# Steady state
ss = hank_1a.hank_ss()

# Compute jacobian of household block
J_ha = jac.get_G([hank_1a.household_trans], ['r', 'w', 'Div', 'Tax'], [], [], T, ss, save=True)

## Using DAG

# define model
exogenous = ['G', 'markup', 'rstar']
unknowns = ['pi', 'w', 'Y']
targets = ['nkpc_res', 'asset_mkt', 'labor_mkt']
outputs = ['pi', 'Tax', 'Y', 'Div', 'r', 'i', 'N', 'w', 'A', 'NS']
block_list=[J_ha, hank_1a.firm, hank_1a.monetary, hank_1a.fiscal, hank_1a.nkpc, hank_1a.mkt_clearing]
dZ = {'Z': 0.9**np.arange(T)}

# run once
curlyJs, required = aux_jac.step0_irf_preliminaries(block_list, dZ, unknowns, ss)
H_U_unpacked = aux_jac.step1_irf_forward_HU(curlyJs, unknowns, targets, required)
J_curlyZ_dZ = aux_jac.step2_irf_forward_Z(curlyJs, dZ, targets, outputs, required)
dU = aux_jac.step3_irf_solve_dU(H_U_unpacked, J_curlyZ_dZ, unknowns, targets, T)
dX = aux_jac.step4_irf_forward_U_combine(curlyJs, dU, J_curlyZ_dZ, outputs, required, T)

# Step 1
res = %timeit -o -q H_U_unpacked = aux_jac.step1_irf_forward_HU(curlyJs, unknowns, targets, required)
times['DAG'][1]['ha'] = np.mean(res.timings)

# Step 2
res = %timeit -o -q J_curlyZ_dZ = aux_jac.step2_irf_forward_Z(curlyJs, dZ, targets, outputs, required)
times['DAG'][2]['ha'] = np.mean(res.timings)

# Step 3
res = %timeit -o -q dU = aux_jac.step3_irf_solve_dU(H_U_unpacked, J_curlyZ_dZ, unknowns, targets, T)
times['DAG'][3]['ha'] = np.mean(res.timings)

# Step 4
res = %timeit -o -q dX = aux_jac.step4_irf_forward_U_combine(curlyJs, dU, J_curlyZ_dZ, outputs, required, T)
times['DAG'][4]['ha'] = np.mean(res.timings)

# Compute total
avg = 0
for i in range(1,5): avg += np.mean(times['DAG'][i]['ha'])
times['DAG'][0]['ha'] = avg

## No DAG

# define model
unknowns = ['Y', 'L', 'w', 'pi', 'Div', 'Tax', 'r']
targets = ['res_L', 'res_Div', 'res_r', 'res_Tax', 'nkpc_res', 'asset_mkt', 'labor_mkt']
block_list=[J_ha, hank_1a.ha1_nodag]

# run once
curlyJs, required = aux_jac.step0_irf_preliminaries(block_list, dZ, unknowns, ss)
H_U_unpacked = aux_jac.step1_irf_forward_HU(curlyJs, unknowns, targets, required)
J_curlyZ_dZ = aux_jac.step2_irf_forward_Z(curlyJs, dZ, targets, outputs, required)
dU = aux_jac.step3_irf_solve_dU(H_U_unpacked, J_curlyZ_dZ, unknowns, targets, T)
dX = aux_jac.step4_irf_forward_U_combine(curlyJs, dU, J_curlyZ_dZ, outputs, required, T)

# Step 1
res = %timeit -o -q H_U_unpacked = aux_jac.step1_irf_forward_HU(curlyJs, unknowns, targets, required)
times['no_DAG'][1]['ha'] = np.mean(res.timings)

# Step 2
res = %timeit -o -q J_curlyZ_dZ = aux_jac.step2_irf_forward_Z(curlyJs, dZ, targets, outputs, required)
times['no_DAG'][2]['ha'] = np.mean(res.timings)

# Step 3
res = %timeit -o -q dU = aux_jac.step3_irf_solve_dU(H_U_unpacked, J_curlyZ_dZ, unknowns, targets, T)
times['no_DAG'][3]['ha'] = np.mean(res.timings)

# Step 4
res = %timeit -o -q dX = aux_jac.step4_irf_forward_U_combine(curlyJs, dU, J_curlyZ_dZ, outputs, required, T)
times['no_DAG'][4]['ha'] = np.mean(res.timings)

# Compute total
avg = 0
for i in range(1,5): avg += np.mean(times['no_DAG'][i]['ha'])
times['no_DAG'][0]['ha'] = avg

Two-asset hank

In [58]:
# Steady state
ss = hank_2a.hank_ss(noisy=False)

# Compute jacobian of household block
J_ha = jac.get_G([hank_2a.household_inc], ['ra', 'rb', 'beta', 'w', 'N', 'tax'], [], [], T, ss)

## Using DAG

# define model
exogenous = ['rstar', 'markup', 'G', 'beta', 'Z', 'rinv_shock', 'markup_w']
unknowns = ['r', 'w', 'Y']
targets = ['asset_mkt', 'fisher', 'wnkpc']
G_pricing = jac.get_G([hank_2a.pricing], ['mc', 'r', 'Y'], ['pi'], ['nkpc'],  T, ss)
G_arbitrage = jac.get_G([hank_2a.arbitrage], ['div', 'r'], ['p'],['equity'], T, ss)
G_production = jac.get_G([hank_2a.labor, hank_2a.investment], ['Y', 'w', 'r', 'Z'], ['Q', 'K'], ['inv', 'val'], T, ss)
block_list = [J_ha, G_pricing, G_arbitrage, G_production, hank_2a.dividend, hank_2a.taylor, hank_2a.fiscal, hank_2a.finance, hank_2a.wage, hank_2a.union, hank_2a.mkt_clearing]
outputs = None
dZ = {'Z': 0.9**np.arange(T)}

# run once
curlyJs, required = aux_jac.step0_irf_preliminaries(block_list, dZ, unknowns, ss)
H_U_unpacked = aux_jac.step1_irf_forward_HU(curlyJs, unknowns, targets, required)
J_curlyZ_dZ = aux_jac.step2_irf_forward_Z(curlyJs, dZ, targets, outputs, required)
dU = aux_jac.step3_irf_solve_dU(H_U_unpacked, J_curlyZ_dZ, unknowns, targets, T)
dX = aux_jac.step4_irf_forward_U_combine(curlyJs, dU, J_curlyZ_dZ, outputs, required, T)

# Step 1
res = %timeit -o -q H_U_unpacked = aux_jac.step1_irf_forward_HU(curlyJs, unknowns, targets, required)
times['DAG'][1]['ha2'] = np.mean(res.timings)

# Step 2
res = %timeit -o -q J_curlyZ_dZ = aux_jac.step2_irf_forward_Z(curlyJs, dZ, targets, outputs, required)
times['DAG'][2]['ha2'] = np.mean(res.timings)

# Step 3
res = %timeit -o -q dU = aux_jac.step3_irf_solve_dU(H_U_unpacked, J_curlyZ_dZ, unknowns, targets, T)
times['DAG'][3]['ha2'] = np.mean(res.timings)

# Step 4
res = %timeit -o -q dX = aux_jac.step4_irf_forward_U_combine(curlyJs, dU, J_curlyZ_dZ, outputs, required, T)
times['DAG'][4]['ha2'] = np.mean(res.timings)

# Compute total
avg = 0
for i in range(1,5): avg += np.mean(times['DAG'][i]['ha2'])
times['DAG'][0]['ha2'] = avg

## No DAG

# define model
unknowns = ['piw', 'psiw', 'rb', 'ra', 'tax', 'i', 'psip', 'I', 'Q', 'w', 'N', 'K', 'div', 'p', 'pi', 'mc', 'r', 'Y']
targets = ['nkpc', 'equity', 'res_N', 'res_mc', 'inv', 'val', 'res_psip', 'res_I', 'res_div', 'res_i', 'res_tax', 'res_rb',
        'res_ra', 'fisher', 'res_piw', 'res_psiw', 'wnkpc', 'asset_mkt']
block_list = [J_ha, hank_2a.ha2_nodag]

# run once
curlyJs, required = aux_jac.step0_irf_preliminaries(block_list, dZ, unknowns, ss)
H_U_unpacked = aux_jac.step1_irf_forward_HU(curlyJs, unknowns, targets, required)
J_curlyZ_dZ = aux_jac.step2_irf_forward_Z(curlyJs, dZ, targets, outputs, required)
dU = aux_jac.step3_irf_solve_dU(H_U_unpacked, J_curlyZ_dZ, unknowns, targets, T)
dX = aux_jac.step4_irf_forward_U_combine(curlyJs, dU, J_curlyZ_dZ, outputs, required, T)

# Step 1
res = %timeit -o -q H_U_unpacked = aux_jac.step1_irf_forward_HU(curlyJs, unknowns, targets, required)
times['no_DAG'][1]['ha2'] = np.mean(res.timings)

# Step 2
res = %timeit -o -q J_curlyZ_dZ = aux_jac.step2_irf_forward_Z(curlyJs, dZ, targets, outputs, required)
times['no_DAG'][2]['ha2'] = np.mean(res.timings)

# Step 3
res = %timeit -o -q dU = aux_jac.step3_irf_solve_dU(H_U_unpacked, J_curlyZ_dZ, unknowns, targets, T)
times['no_DAG'][3]['ha2'] = np.mean(res.timings)

# Step 4
res = %timeit -o -q dX = aux_jac.step4_irf_forward_U_combine(curlyJs, dU, J_curlyZ_dZ, outputs, required, T)
times['no_DAG'][4]['ha2'] = np.mean(res.timings)

# Compute total
avg = 0
for i in range(1,5): avg += np.mean(times['no_DAG'][i]['ha2'])
times['no_DAG'][0]['ha2'] = avg
In [59]:
if savedata: np.save('import_export/results/irf_speed',[times])

7.4 Nonlinear impulse responses

In [60]:
times = {}
T = 300
dZ = 0.01*0.8**np.arange(T)

# Krusell-Smith
ss = ks.ks_ss()
unknowns=['K']
targets=['asset_mkt']
block_list=[ks.firm, ks.mkt_clearing, ks.household]
G = jac.get_G(block_list=block_list, exogenous=unknowns, unknowns=[], targets=[], outputs=targets, T=T, ss=ss)
H_U = jac.pack_jacobians(G, inputs=unknowns, outputs=targets, T=T)
td_nonlinear = nonlinear.td_solve(ss, block_list, unknowns, targets, H_U=H_U, Z=ss['Z']*(1+dZ), noisy=False)
res = %timeit -o -q td_nonlinear = nonlinear.td_solve(ss, block_list, unknowns, targets, H_U=H_U, Z=ss['Z']*(1+dZ), noisy=False)
times['ks'] = np.mean(res.timings)

# Krusell-Smith HD
ss = ks.ks_ss(nS=50, nA=5000)
unknowns=['K']
targets=['asset_mkt']
block_list=[ks.firm, ks.mkt_clearing, ks.household]
G = jac.get_G(block_list=block_list, exogenous=unknowns, unknowns=[], targets=[], outputs=targets, T=T, ss=ss)
H_U = jac.pack_jacobians(G, inputs=unknowns, outputs=targets, T=T)
td_nonlinear = nonlinear.td_solve(ss, block_list, unknowns, targets, H_U=H_U, Z= ss['Z']*(1+dZ), noisy=False)
res = %timeit -o -q td_nonlinear = nonlinear.td_solve(ss, block_list, unknowns, targets, H_U=H_U, Z=ss['Z']*(1+dZ), noisy=False)
times['hd_ks'] = np.mean(res.timings)

# 1-asset hank
ss = hank_1a.hank_ss()
unknowns = ['pi', 'w', 'Y']
targets = ['nkpc_res', 'asset_mkt', 'labor_mkt']
block_list = [hank_1a.household_trans, hank_1a.firm, hank_1a.monetary, hank_1a.fiscal, hank_1a.nkpc, hank_1a.mkt_clearing]
G = jac.get_G(block_list=block_list, exogenous=unknowns, unknowns=[], targets=[], outputs=targets, T=T, ss=ss)
H_U = jac.pack_jacobians(G, inputs=unknowns, outputs=targets, T=T)
td_nonlinear = nonlinear.td_solve(ss, block_list, unknowns, targets, H_U=H_U, Z= ss['Z']*(1+dZ), noisy=False)
res = %timeit -o -q td_nonlinear = nonlinear.td_solve(ss, block_list, unknowns, targets, H_U=H_U, Z=ss['Z']*(1+dZ), noisy=False)
times['ha'] = np.mean(res.timings)

# 2-asset hank
ss = hank_2a.hank_ss(noisy=False)
unknowns = ['r', 'w', 'Y', 'K', 'Q', 'pi', 'p']
targets = ['asset_mkt', 'fisher', 'wnkpc', 'inv', 'val', 'equity', 'nkpc']
block_list = [hank_2a.household_inc, hank_2a.pricing, hank_2a.arbitrage, hank_2a.labor, hank_2a.investment, hank_2a.dividend, hank_2a.taylor, hank_2a.fiscal, hank_2a.finance, hank_2a.wage, hank_2a.union, hank_2a.mkt_clearing]
G = jac.get_G(block_list=block_list, exogenous=unknowns, unknowns=[], targets=[], outputs=targets, T=T, ss=ss)
H_U = jac.pack_jacobians(G, inputs=unknowns, outputs=targets, T=T)
td_nonlinear = nonlinear.td_solve(ss, block_list, unknowns, targets, H_U=H_U, Z=ss['Z']*(1+dZ), noisy=False)
res = %timeit -o -q td_nonlinear = nonlinear.td_solve(ss, block_list, unknowns, targets, H_U=H_U, Z=ss['Z']*(1+dZ), noisy=False)
times['ha2'] = np.mean(res.timings)

if savedata: np.save('import_export/results/nonlinear_speed', [times])