[Python] IRFs Appear Reversed in RBC Model Simulation - Problem with Recursive Shock...

I'm trying to simulate a fairly standard RBC model. The equations used in the code are all correct based on the chosen specification, so my question relates more to the implementation side.

Specifically, I want to use the technology shock generated in z_hat:

# Generate the autoregressive process for the technology shock:
for t in range(1, T):
# Each period the shock decays by the factor rho_z, mimicking an AR(1) process.
z_hat[t] = rho_z * z_hat[t - 1]


to influence the other model equations recursively. The responses of the model variables are captured in the for loop that follows.

The issue is this: when I run the code as it is, the impulse response functions (IRFs) appear to be flipped. I only get the expected graph if I reverse the order of the calculated IRFs for each variable.

Has anyone encountered a similar issue or can spot what might be going wrong in how the IRFs are being generated or stored? I attached the code which created the IRFs and the full code.

Thanks in advance!

Important part of the code:

# Generate the autoregressive process for the technology shock:
for t in range(1, T):
# Each period the shock decays by the factor rho_z, mimicking an AR(1) process.
z_hat[t] = rho_z * z_hat[t - 1]

# Initialize arrays to store the log-deviations from steady state for various variables.
# The deviations represent percentage deviations from the benchmark steady state.
c_hat = np.zeros(T+1) # Consumption log-deviation, extra period for terminal condition.
i_hat = np.zeros(T) # Investment log-deviation.
k_hat = np.zeros(T+1) # Capital log-deviation, extra period due to accumulation.
y_hat = np.zeros(T) # Output log-deviation.
r_hat = np.zeros(T) # Interest rate log-deviation.
w_hat = np.zeros(T) # Wage log-deviation.
l_hat = np.zeros(T) # Labor log-deviation.

# Terminal condition is assumed for consumption: c_T = 0 (i.e., no deviation at the final period).

# Set initial conditions: assume no deviation in capital and consumption at time t = 0.
k_hat[0] = 0
c_hat[0] = 0

# Loop over each time period to simulate the dynamic responses (Impulse Response Functions - IRFs)
for t in range(T):

# 1. Factor Prices from Marginal Products:
# The wage (w_hat) is determined as the sum of the technology shock and contributions from
# capital and labor deviations scaled by the capital share.
w_hat[t] = z_hat[t] + alpha * k_hat[t] - alpha * l_hat[t]
# The rental rate (r_hat) is calculated using the marginal product conditions;
# it is affected by the deviation of capital and labor.
r_hat[t] = (alpha - 1) * k_hat[t] + z_hat[t] + (1 - alpha) * l_hat[t]

# 2. Labor Supply from the Intratemporal Condition:
# The log-deviation of labor (l_hat) is updated based on the difference between wages and consumption,
# scaled by the ratio ((1-L_ss)/L_ss), which captures the elasticity in this specification.
l_hat[t] = (w_hat[t] - c_hat[t]) * ((1 - L_ss) / L_ss)

# 3. Output Computation:
# Output is determined by combining the contributions of capital and labor deviations weighted by alpha and
# (1-alpha), respectively, then adding the technology shock.
# This update is performed after updating the labor deviation.
y_hat[t] = alpha * k_hat[t] + (1 - alpha) * l_hat[t] + z_hat[t]

# 4. Investment from the Resource Constraint:
# Investment log-deviation is calculated using the resource constraint where output minus the scaled consumption
# is divided by the share of investment in output.
i_hat[t] = (y_hat[t] - phi_c * c_hat[t]) / phi_i

# 5. Capital Accumulation:
# The next period's capital log-deviation is determined by the current capital after depreciation and the
# new investment weighted by delta, representing the adjustment process.
k_hat[t + 1] = (1 - delta) * k_hat[t] + delta * i_hat[t]

# 6. Euler Equation (Consumption Dynamics):
# According to the Euler equation in a log-linearized setting, the change in the consumption deviation
# equals the product of the steady state return factor (beta*q_ss) and the rental rate deviation.
c_hat[t + 1] = c_hat[t] + (beta * q_ss) * r_hat[t]


Full code:

###############################################################################
# Loading Required Packages
###############################################################################
# Import the numpy package for efficient numerical computing and array manipulation.
import numpy as np
# Import sympy for symbolic mathematics (e.g., symbolic algebra, calculus).
import sympy as sp
# Import matplotlib's pyplot for plotting and creating visualizations.
import matplotlib.pyplot as plt
# Set the seed for NumPy's random number generator to ensure that any random processes produce
# the same output every time the code is run (this aids in reproducibility).
np.random.seed(42) # optional, for reproducibility
# Import the math module for accessing standard mathematical functions.
import math

###############################################################################
# End of package loading
###############################################################################


###############################################################################
###############################################################################
# 1% Shock in Technology: Set up the DSGE model environment and steady state.
###############################################################################
###############################################################################

# Starting parameter values for our DSGE model:

# beta: Discount factor, reflecting how future payoffs are weighted relative to present ones.
beta = 0.96
# delta: Depreciation rate of capital per period.
delta = 0.1
# p: Capital's share of output (alternatively represented as alpha in many models).
p = 0.36
# Lss: Steady-state (or initial) labor input; here, representing a small fraction for the population.
Lss = 0.01 # Population starting value
# z: Technology level; initially set to a small value (here equivalent to a 1% initial shock).
z = 0.01

# Define functions to compute steady state values using the model's equations.
# These functions calculate key economic quantities as functions of the parameters and z.

def K_fun(beta, delta, p, Lss, z):
"""
Computes the steady state level of Capital (K) using the model's formula.
The expression is derived from equating marginal products and other steady state conditions.
"""
global K_expr
# The expression involves the technology term z, labor Lss and other parameters.
K_expr = (p * Lss**(1 - p) * z / (1/beta - 1 + delta))**(1/(1 - p))
return K_expr

def q_fun(beta, delta, p, Lss, z):
"""
Computes the steady state marginal product of capital (q) which is related
to the return on capital.
"""
global q_expr
# The function inverts part of the K_fun expression and rescales it by z, p and Lss.
q_expr = z * p * (p * Lss**(1 - p) * z / (1/beta - 1 + delta))**(-1) * Lss**(1 - p)
return q_expr

def w_fun(beta, delta, p, Lss, z):
"""
Computes the steady state wage (w) as given by the marginal product of labor.
"""
global w_expr
# Note the use of the exponent p/(1-p) in the wage formula, capturing the relative power
# of capital’s effect adjusted by the labor share.
w_expr = z * (1 - p) * Lss**(-p) * (p * Lss**(1 - p) * z / (1/beta - 1 + delta))**(p/(1 - p))
return w_expr

def Y_fun(beta, delta, p, Lss, z):
"""
Computes the steady state level of Output (Y) using the production function.
This function assumes a Cobb-Douglas production function depending on K and L.
"""
global Y_expr
# K_expr is used here, and it is assumed to be computed before using K_fun.
Y_expr = z * (K_expr)**p * Lss**(1 - p)
return Y_expr

def I_fun(beta, delta, p, Lss, z):
"""
Computes the level of Investment (I) in the steady state.
Investment is proportional to the capital stock through the depreciation rate.
"""
global I_expr
I_expr = delta * K_expr
return I_expr

def C_fun(beta, delta, p, Lss, z):
"""
Computes steady state Consumption (C) as the residual of output after investment.
"""
global C_expr
C_expr = Y_expr - I_expr
return C_expr

def L_fun(beta, delta, p, Lss, z):
"""
Provides the value of steady state labor. Here it is set to L_ss (previously defined),
implying fixed labor input in the model.
"""
global L_expr
L_expr = Lss
return L_expr

# ---------------------------
# 2. Define the steady state equations (as expressions) in terms of z
# ---------------------------
# Compute the steady state values using the defined functions.
# These are the benchmark values around which log-linearized deviations are computed.
K_val = K_fun(beta, delta, p, Lss, z)
q_val = q_fun(beta, delta, p, Lss, z)
w_val = w_fun(beta, delta, p, Lss, z)
Y_val = Y_fun(beta, delta, p, Lss, z)
I_val = I_fun(beta, delta, p, Lss, z)
C_val = C_fun(beta, delta, p, Lss, z)
L_val = L_fun(beta, delta, p, Lss, z)

# The above assignments are immediately repeated to store the steady state as "_ss" variables,
# which might be used later to make the notation clearer that these are steady-state benchmarks.
K_ss = K_fun(beta, delta, p, Lss, z)
q_ss = q_fun(beta, delta, p, Lss, z)
w_ss = w_fun(beta, delta, p, Lss, z)
Y_ss = Y_fun(beta, delta, p, Lss, z)
I_ss = I_fun(beta, delta, p, Lss, z)
C_ss = C_fun(beta, delta, p, Lss, z)
L_ss = L_fun(beta, delta, p, Lss, z)

# Print the steady state values to verify the baseline of the model. These are calculated at z=1.
print("Steady State Values (at z=1):")
print("K_ss =", K_val)
print("q_ss =", q_val)
print("w_ss =", w_val)
print("Y_ss =", Y_val)
print("I_ss =", I_val)
print("C_ss =", C_val)
print("L_ss =", L_val)

# Define additional parameters for the simulation:
alpha = p # alpha is set equal to capital share, which in some literature uses alpha.
T = 20 # T defines the total number of time periods for the IRF simulation.
rho_z = 0.95 # Autoregressive parameter for the technology shock process.

# From the steady state, compute additional model parameters:
# r_ss is the steady state real interest rate derived from the Euler equation.
r_ss = 1 / beta - 1 + delta
# phi_c and phi_i are ratios of consumption and investment to output in the steady state,
# and are used to scale the IRF responses.
phi_c = C_ss / Y_ss
phi_i = I_ss / Y_ss

# Initialize the technology shock process:
# z_hat is an array representing the shock process over T periods.
z_hat = np.zeros(T)
z_hat[0] = 1 # At time 0, we introduce a 1% technology shock.

# Generate the autoregressive process for the technology shock:
for t in range(1, T):
# Each period the shock decays by the factor rho_z, mimicking an AR(1) process.
z_hat[t] = rho_z * z_hat[t - 1]

# Initialize arrays to store the log-deviations from steady state for various variables.
# The deviations represent percentage deviations from the benchmark steady state.
c_hat = np.zeros(T+1) # Consumption log-deviation, extra period for terminal condition.
i_hat = np.zeros(T) # Investment log-deviation.
k_hat = np.zeros(T+1) # Capital log-deviation, extra period due to accumulation.
y_hat = np.zeros(T) # Output log-deviation.
r_hat = np.zeros(T) # Interest rate log-deviation.
w_hat = np.zeros(T) # Wage log-deviation.
l_hat = np.zeros(T) # Labor log-deviation.

# Terminal condition is assumed for consumption: c_T = 0 (i.e., no deviation at the final period).

# Set initial conditions: assume no deviation in capital and consumption at time t = 0.
k_hat[0] = 0
c_hat[0] = 0

# Loop over each time period to simulate the dynamic responses (Impulse Response Functions - IRFs)
for t in range(T):

# 1. Factor Prices from Marginal Products:
# The wage (w_hat) is determined as the sum of the technology shock and contributions from
# capital and labor deviations scaled by the capital share.
w_hat[t] = z_hat[t] + alpha * k_hat[t] - alpha * l_hat[t]
# The rental rate (r_hat) is calculated using the marginal product conditions;
# it is affected by the deviation of capital and labor.
r_hat[t] = (alpha - 1) * k_hat[t] + z_hat[t] + (1 - alpha) * l_hat[t]

# 2. Labor Supply from the Intratemporal Condition:
# The log-deviation of labor (l_hat) is updated based on the difference between wages and consumption,
# scaled by the ratio ((1-L_ss)/L_ss), which captures the elasticity in this specification.
l_hat[t] = (w_hat[t] - c_hat[t]) * ((1 - L_ss) / L_ss)

# 3. Output Computation:
# Output is determined by combining the contributions of capital and labor deviations weighted by alpha and
# (1-alpha), respectively, then adding the technology shock.
# This update is performed after updating the labor deviation.
y_hat[t] = alpha * k_hat[t] + (1 - alpha) * l_hat[t] + z_hat[t]

# 4. Investment from the Resource Constraint:
# Investment log-deviation is calculated using the resource constraint where output minus the scaled consumption
# is divided by the share of investment in output.
i_hat[t] = (y_hat[t] - phi_c * c_hat[t]) / phi_i

# 5. Capital Accumulation:
# The next period's capital log-deviation is determined by the current capital after depreciation and the
# new investment weighted by delta, representing the adjustment process.
k_hat[t + 1] = (1 - delta) * k_hat[t] + delta * i_hat[t]

# 6. Euler Equation (Consumption Dynamics):
# According to the Euler equation in a log-linearized setting, the change in the consumption deviation
# equals the product of the steady state return factor (beta*q_ss) and the rental rate deviation.
c_hat[t + 1] = c_hat[t] + (beta * q_ss) * r_hat[t]

# The consumption and capital paths are reversed to align the time axis for plotting,
# recognizing that c(t+1) and k(t+1) are plotted rather than c(t) and k(t).
c_hat_aligned = np.insert(c_hat[::-1][:-1], 0, 0.0) # Insert a zero at the beginning after reversing.
k_hat_aligned = np.insert(k_hat[::-1][:-1], 0, 0.0) # Similarly for capital.

# Prepare the variables for plotting by converting log-deviations into percentage deviations.
# The dictionary "variables" maps variable names to their corresponding reversed IRF paths.
variables = {
'Output (Y)': y_hat[::-1],
'Capital (K)': k_hat_aligned,
'Investment (I)': i_hat[::-1],
'Consumption (C)': c_hat_aligned,
'Wage (w)': w_hat[::-1],
'Interest Rate (r)': r_hat[::-1],
'Labor (L)': l_hat[::-1],
}

# Set up the plot size for the IRFs.
plt.figure(figsize=(12, 6))

# Plot each variable's deviation over time.
for name, var in variables.items():
# Create a time array that matches the length of the variable series.
time_var = np.arange(len(var))
plt.plot(time_var, var, label=name)

# Also plot the technology shock (z_hat) for reference.
plt.plot(np.arange(len(z_hat)), z_hat, label='Technology (z)', linewidth=1, color='black')

# Add a horizontal reference line at zero to clearly see deviations from the steady state.
plt.axhline(0, color='gray', linestyle='--')
# Add plot title and axis labels.
plt.title("IRFs from Log-Linearized DSGE Model (Shock at t = 0)")
plt.xlabel("Time")
plt.ylabel("% Deviation from Steady State")
# Add a legend to distinguish between different variables.
plt.legend()
# Display grid lines for improved readability.
plt.grid(True)
# Adjust the layout to ensure nothing is clipped.
plt.tight_layout()
# Render the final plot.
plt.show()

Continue reading...