In [1]:
import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt

import scipy.linalg

sns.set_style("whitegrid")


#### Coefficients¶

In [2]:
alpha = 1 / 3
beta = .99
delta = .025
psi = 2.
rho = .98
h = 1 / 3

# in following variable names: 'ss ~ steady state'


### Compute Objects¶

In [3]:
ss_capital = h * ((1/beta - 1 + delta) / alpha) ** (1 / (alpha - 1))
ss_prod = (ss_capital ** alpha) * (h ** (1 - alpha))
ss_consum = ss_prod - delta * ss_capital
ss_invest = delta * ss_capital

In [4]:
m = 1 - beta * (1 - delta)
b = alpha + (1-delta) * (ss_capital / ss_prod)

A = np.array([
[m*(1-alpha)         , -m     , 1, -m*(1-alpha)],
[0                   , 0      , 0, 0           ],
[ss_capital / ss_prod, 0      , 0, 0           ],
[0                   , 1 / rho, 0, 0           ],
])

B = np.array([
[0    , 0, 1                   , 0             ],
[alpha, 1, -1                  , -(psi + alpha)],
[b    , 1, -ss_consum / ss_prod, 1 - alpha     ],
[0    , 1, 0                   , 0             ],
])


### g)¶

In the function scipy.linalg.ordqz I set sort="ouc" which stands for outside unit circle. This achieves the desired order for our problem.

In [5]:
S, T, _, _, Q, Z =  scipy.linalg.ordqz(A, B, sort="ouc")


### h)¶

In [6]:
eigen_values = np.diag(T) / np.diag(S)

I = np.where(np.abs(eigen_values) < 1)[0]

ns = len(I)

S_11, T_11, Z_11 = (m[:ns, :ns] for m in (S, T, Z))
Z_21 = Z[ns:, :ns]

Z_11_inverse = np.linalg.inv(Z_11)


#### (Linear) Mappings from States to States and States to Controls:¶

In [7]:
S2S = Z_11 @ np.linalg.inv(S_11) @ T_11 @ np.linalg.inv(Z_11) # states to states
S2C = Z_21 @ np.linalg.inv(Z_11)                              # states to controls


### j)¶

In [8]:
def simulate_impulse_response(S2S, S2C, shock_level, n_periods):
"""Simulate impulse response path.

Args:
S2S (np.ndarray): Array mapping states to states
S2C (np.ndarray): Array mapping states to controls
shock_level (float): Size of technology shock
n_periods (int): Number of periods to simulate.

Returns:
df (pd.DataFrame): Data frame with simulated variables.
Columns are 'periods', 'variable' and 'value', where
'variable' is in "capital", "productivity",
"consumption", "labor", "investment", "prodution"

"""
ns = len(S2S)
nc = len(S2C)

# result will contain state and controls for all time periods plus further info
# columns of result will be [capital, productivity, consum, labor, investment,
# production]
result = np.empty((n_periods, ns + nc + 2))

state = np.array([0, shock_level])

for t in range(n_periods):

controls = S2C @ state
updated_state = S2S @ state

result[t, :ns] = state
result[t, ns:(ns + nc)] = controls

investment = (updated_state[0] - (1 - delta) * state[0]) / delta
prod = ss_consum / ss_prod * controls[0] + ss_invest / ss_prod * investment

result[t, -2] = investment
result[t, -1] = prod

state = updated_state

result /= shock_level

df = pd.DataFrame(result, columns=[
"capital", "productivity", "consumption", "labor", "investment", "prodution"
])
df.index.name = "periods"
df = df.reset_index().melt(id_vars="periods", var_name="variable")
return df


#### Actual Computation:¶

In [9]:
impulse_response = simulate_impulse_response(
S2S, S2C, shock_level=.007, n_periods=400
)

In [10]:
fig, ax = plt.subplots(1)
fig.set_size_inches(15, 8)

sns.lineplot(
x="periods",
y="value",
data=impulse_response,
hue="variable",
legend="brief",
linewidth=2,
ax=ax,
)

ax.set(xlabel="Periods")
ax.set_frame_on(False)
ax.xaxis.label.set_size(13)
ax.yaxis.label.set_size(13)
plt.legend(prop={'size': 13})
plt.show()