```python
import numpy as np
import pandas as pd
import linearsolve as ls
import matplotlib.pyplot as plt
plt.style.use('classic')
%matplotlib inline
```

# Class 18: Exogenous Shocks in the New-Keynesian Model

In this notebook, we will compute impulse responses to exogenous shocks in the new-Keynesian model and interpret the results. Then we compute a stochastic simulation of the model.

## The New-Keynesian Model

The most basic version of the New-Keynesian Model can be expressed as:

\begin{align}
y_t & = E_t y_{t+1}  - \left( r_{t} - \bar{r}\right) + g_t\\
i_{t} & = r_{t} + E_t \pi_{t+1}\\
i_{t} & = \bar{r} + \pi^T + \phi_{\pi}\big(\pi_t - \pi^T\big) + \phi_{y}\big(y_t - \bar{y}\big) + v_t\\
\pi_t -\pi^T & = \beta \left( E_t\pi_{t+1} - \pi^T\right)  + \kappa (y_t -\bar{y})+ u_t,
\end{align}

where: $y_t$ is (log) output, $r_t$ is the real interest rate, $i_t$ is the nominal interest rate, $\pi_t$ is the rate of inflation between periods $t-1$ and $t$, $\bar{r}$ is the long-run average real interest rate or the *natural rate of interest*, $\beta$ is the household's subjective discount factor, and $\pi^T$ is the central bank's inflation target. The coeffieints $\phi_{\pi}$ and $\phi_{y}$ reflect the degree of intensity to which the central bank *endogenously* adjusts the nominal interest rate in response to movements in inflation and output.

The variables $g_t$, $u_t$, and $v_t$ represent exogenous shocks to aggregate demand, inflation, and monetary policy. They follow AR(1) processes:

\begin{align}
g_{t+1} & = \rho_g g_{t} + \epsilon^g_{t+1}\\
u_{t+1} & = \rho_u u_{t} + \epsilon^u_{t+1}\\
v_{t+1} & = \rho_v v_{t} + \epsilon^v_{t+1}.
\end{align}

The goal is to compute impulse responses in the model to all three exogenous shocks. We will use the following parameterization:

| $$\bar{y}$$ | $$\beta$$ | $$\bar{r}$$    | $$\kappa$$ | $$\pi^T$$ | $$\phi_{\pi}$$ | $$\phi_y$$ | $$\rho_g$$ | $$\rho_u$$ | $$\rho_v$$ | 
|-------------|-----------|----------------|------------|-----------|----------------|------------|------------|------------|------------|
| 0           | 0.995     | $$-\log\beta$$ | 0.25       | 0.02/4    | 1.5            | 0.5/4      | 0.5        | 0.5        | 0.5        |

## Solution


```python
# Create a variable called 'parameters' that stores the model parameter values in a Pandas Series. CELL PROVIDED
parameters = pd.Series(dtype=float)
parameters['y_bar'] = 0
parameters['beta'] = 0.995
parameters['r_bar'] = -np.log(parameters.beta)
parameters['kappa'] = 0.25
parameters['pi_T'] = 0.02/4
parameters['phi_pi'] = 1.5
parameters['phi_y'] = 0.5/4
parameters['rho_g'] = 0.5
parameters['rho_u'] = 0.5
parameters['rho_v'] = 0.5

# Create variable called 'var_names' that stores the variable names in a list with state variables ordered first
var_names = ['g','u','v','y','pi','i','r']

# Create variable called 'shock_names' that stores an exogenous shock name for each state variable.
shock_names = ['e_g','e_u','e_v']

# Define a function that evaluates the equilibrium conditions of the model solved for zero.
def equilibrium_equations(variables_forward,variables_current,parameters):
    
    # Parameters
    p = parameters
    
    # Current variables
    cur = variables_current
    
    # Forward variables
    fwd = variables_forward

    # IS equation
    is_equation = fwd.y - (cur.r -p.r_bar) + cur.g - cur.y
    
    # Fisher_equation
    fisher_equation = cur.r + fwd.pi - cur.i
    
    # Monetary policy
    monetary_policy = p.r_bar + p.pi_T + p.phi_pi*(cur.pi - p.pi_T) + p.phi_y*cur.y + cur.v - cur.i
    
    # Phillips curve
    phillips_curve = p.beta*(fwd.pi- p.pi_T) + p.kappa*cur.y + cur.u - (cur.pi-p.pi_T)
    
    # Demand process
    demand_process = p.rho_g*cur.g - fwd.g
    
    # Monetary policy process
    monetary_policy_process = p.rho_v*cur.v - fwd.v
    
    # Inflation process
    inflation_process = p.rho_u*cur.u - fwd.u
    
    
    # Stack equilibrium conditions into a numpy array
    return np.array([
            is_equation,
            fisher_equation,
            monetary_policy,
            phillips_curve,
            demand_process,
            monetary_policy_process,
            inflation_process
        ])

# Initialize the model into a variable named 'nk_model'
nk_model = ls.model(equations = equilibrium_equations,
                 n_states=3,
                 var_names=var_names,
                 shock_names=shock_names,
                 parameters = parameters)


# Compute the steady state numerically using .compute_ss() method of nk_model
guess = [0,0,0,0,0.01,0.01,0.01]
nk_model.compute_ss(guess)
```


```python
# Find the log-linear approximation around the non-stochastic steady state and solve using .approximate_and_solve() method of nk_model
# set argumement 'log_linear' to False
nk_model.approximate_and_solve(log_linear=False)

# Print the solution to the model
print(nk_model.solved())
```

    Solution to the linear system:
    
    g[t+1] = 0.5·g[t]+e_g[t+1]
    
    u[t+1] = 0.5·u[t]+e_u[t+1]
    
    v[t+1] = 0.5·v[t]+e_v[t+1]
    
      y[t] = 0.8909·g[t]-1.7729·u[t]-0.8909·v[t]
    
     pi[t] = 0.4432·g[t]+1.108·u[t]-0.4432·v[t]
    
      i[t] = 0.7762·g[t]+1.4404·u[t]+0.2238·v[t]
    
      r[t] = 0.5546·g[t]+0.8864·u[t]+0.4454·v[t]


**Questions** 

1. Will a positive shock to $g_t$ increase or decrease $y_t$?
2. Will a positive shock to $g_t$ increase or decrease $\pi_t$?
3. Will a positive shock to $u_t$ increase or decrease $y_t$?
4. Will a positive shock to $u_t$ increase or decrease $\pi_t$?

**Answers**

1. The coefficient on $g_t$ in the solution for $y_t$ is *positive* so the shock will *increase* $y_t$ <!-- answer -->
2. The coefficient on $g_t$ in the solution for $\pi_t$ is *positive* so the shock will *increase* $\pi_t$ <!-- answer -->
1. The coefficient on $u_t$ in the solution for $y_t$ is *negative* so the shock will *decrease* $y_t$ <!-- answer -->
4. The coefficient on $u_t$ in the solution for $\pi_t$ is *positive* so the shock will *increase* $\pi_t$ <!-- answer -->

## Impulse Responses

Compute a 21 period impulse response of the model's variables to a:

* 0.01 unit shock to $g_t$
* 0.01/4 unit shock to $u_t$
* 0.01/4 unit shock to $v_t$


```python
# Compute impulse responses
nk_model.impulse(T=21,t0=5,shocks=[0.01,0.01/4,0.01/4])
```

### Demand Shock

Plot the computed impulses responses of the nominal interest rate, the real interest rate, output, and inflation to the 0.01 unit demand shock. Express inflation and interest rates in *annualized* (e.g., multiplied by 4) terms.


```python
# Create figure. CELL PROVIDED
fig = plt.figure(figsize=(12,8))

# Create upper-left axis
ax1 = fig.add_subplot(2,2,1)
# Create upper-right axis
ax2 = fig.add_subplot(2,2,2)
# Create lower-left axis
ax3 = fig.add_subplot(2,2,3)
# Create lower-right axis
ax4 = fig.add_subplot(2,2,4)

# Set axis 1 ylabel
ax1.set_ylabel('% dev from steady state')
# Set axis 2 ylabel
ax2.set_ylabel('% dev from steady state')
# Set axis 3 ylabel
ax3.set_ylabel('% dev from steady state')
# Set axis 4 ylabel
ax4.set_ylabel('% dev from steady state')

# Set axis 1 limits 
ax1.set_ylim([-0.5,3.5])
# Set axis 2 limits
ax2.set_ylim([-0.5,3.5])
# Set axis 3 limits
ax3.set_ylim([-0.5,2])
# Set axis 4 limits
ax4.set_ylim([-0.5,2])

# Plot the nominal interest rate, real interest rate, output, and inflation
(nk_model.irs['e_g']['i']*400).plot(ax=ax1,lw=4,alpha=0.75,title='Nominal Interest (annualized)',grid=True)
(nk_model.irs['e_g']['r']*400).plot(ax=ax2,lw=4,alpha=0.75,title='Real Interest (annualized)',grid=True)
(nk_model.irs['e_g']['y']*100).plot(ax=ax3,lw=4,alpha=0.75,title='Output',grid=True)
(nk_model.irs['e_g']['pi']*400).plot(ax=ax4,lw=4,alpha=0.75,title='Inflation (annualized)',grid=True)
```

### Inflation Shock

Plot the computed impulses responses of the nominal interest rate, the real interest rate, output, and inflation to the 0.01/4 unit inflation shock. Express inflation and interest rates in *annualized* (e.g., multiplied by 4) terms.


```python
# Create figure. CELL PROVIDED
fig = plt.figure(figsize=(12,8))

# Create upper-left axis
ax1 = fig.add_subplot(2,2,1)
# Create upper-right axis
ax2 = fig.add_subplot(2,2,2)
# Create lower-left axis
ax3 = fig.add_subplot(2,2,3)
# Create lower-right axis
ax4 = fig.add_subplot(2,2,4)

# Set axis 1 ylabel
ax1.set_ylabel('% dev from steady state')
# Set axis 2 ylabel
ax2.set_ylabel('% dev from steady state')
# Set axis 3 ylabel
ax3.set_ylabel('% dev from steady state')
# Set axis 4 ylabel
ax4.set_ylabel('% dev from steady state')

# Set axis 1 limits 
ax1.set_ylim([-.5,2])
# Set axis 2 limits
ax2.set_ylim([-0.5,2])
# Set axis 3 limits
ax3.set_ylim([-0.5,0.5])
# Set axis 4 limits
ax4.set_ylim([-0.5,1.5])

# Plot the nominal interest rate, real interest rate, output, and inflation
(nk_model.irs['e_u']['i']*400).plot(ax=ax1,lw=4,alpha=0.75,title='Nominal Interest (annualized)',grid=True)
(nk_model.irs['e_u']['r']*400).plot(ax=ax2,lw=4,alpha=0.75,title='Real Interest (annualized)',grid=True)
(nk_model.irs['e_u']['y']*100).plot(ax=ax3,lw=4,alpha=0.75,title='Output',grid=True)
(nk_model.irs['e_u']['pi']*400).plot(ax=ax4,lw=4,alpha=0.75,title='Inflation (annualized)',grid=True)
```

### Monetary Policy Shock

Plot the computed impulses responses of the nominal interest rate, the real interest rate, output, and inflation to the 0.01/4 unit monetary policy shock. Express inflation and interest rates in *annualized* (e.g., multiplied by 4) terms.


```python
# Create figure. CELL PROVIDED
fig = plt.figure(figsize=(12,8))

# Create upper-left axis
ax1 = fig.add_subplot(2,2,1)
# Create upper-right axis
ax2 = fig.add_subplot(2,2,2)
# Create lower-left axis
ax3 = fig.add_subplot(2,2,3)
# Create lower-right axis
ax4 = fig.add_subplot(2,2,4)

# Set axis 1 ylabel
ax1.set_ylabel('% dev from steady state')
# Set axis 2 ylabel
ax2.set_ylabel('% dev from steady state')
# Set axis 3 ylabel
ax3.set_ylabel('% dev from steady state')
# Set axis 4 ylabel
ax4.set_ylabel('% dev from steady state')

# Set axis 1 limits 
ax1.set_ylim([-0.2,0.6])
# Set axis 2 limits
ax2.set_ylim([-0.2,0.6])
# Set axis 3 limits
ax3.set_ylim([-0.5,0.1])
# Set axis 4 limits
ax4.set_ylim([-0.5,0.1])

# Plot the nominal interest rate, real interest rate, output, and inflation
(nk_model.irs['e_v']['i']*400).plot(ax=ax1,lw=4,alpha=0.75,title='Nominal Interest (annualized)',grid=True)
(nk_model.irs['e_v']['r']*400).plot(ax=ax2,lw=4,alpha=0.75,title='Real Interest (annualized)',grid=True)
(nk_model.irs['e_v']['y']*100).plot(ax=ax3,lw=4,alpha=0.75,title='Output',grid=True)
(nk_model.irs['e_v']['pi']*400).plot(ax=ax4,lw=4,alpha=0.75,title='Inflation (annualized)',grid=True)
```

## Example: Stochastic Simulation

Compute a 401 period stochastic simulation of the model.  

| $$\rho_g$$ | $$\sigma_g$$ | $$\rho_u$$ | $$\sigma_u$$ | $$\rho_v$$ | $$\sigma_v$$ |
|------------|--------------|------------|--------------|------------|--------------|
| 0.5        | 0.015        | 0.5        | 0.003        | 0.5        | 0.001        |

The parameter values are chosen so that the following simulated values are on par with the historical data for the US:

* standard deviation of output
* standard deviation of inflation
* standard deviation of the nominal interest rate
* correlation between output and inflation
* correlation between output and the nominal interest rate


```python
# Form the covariance matrix for the exogenous white noise processes 
Sigma = [[0.015**2,0,0],[0,0.003**2,0],[0,0,0.001**2]]
```


```python
# Compute the stochastic simulation using the .stoch_sim() method of nk_model. Set seed to 126.
nk_model.stoch_sim(T=201,cov_mat=Sigma,seed=126)
```


```python
# Create 12x4 figure. PROVIDED
fig = plt.figure(figsize=(12,12))

# Create axis. PROVIDED
ax1 = fig.add_subplot(3,1,1)

# Plot the output gap
ax1.plot(nk_model.simulated['y']*100)
ax1.set_title('Output gap')
ax1.set_ylabel('Percent')
ax1.grid()

# Create axis. PROVIDED
ax2 = fig.add_subplot(3,1,2)

# Plot inflation (annualized)
ax2.plot(nk_model.simulated['pi']*400)
ax2.set_title('Inflation (annualized)')
ax2.set_ylabel('Percent')
ax2.grid()

# Create axis. PROVIDED
ax3 = fig.add_subplot(3,1,3)

# Plot the nominal interest rate (annualized)
ax3.plot(nk_model.simulated['i']*400)
ax3.set_title('Nominal interest rate (annualized)')
ax3.set_ylabel('Percent')
ax3.grid()
```


```python
# Compute standard deviations of output (times 100), inflation (times 400), and the nominal interest rate (times 400)
nk_model.simulated[['y','pi','i']].std()*[100,400,400]
```




    y     1.686375
    pi    3.779509
    i     6.242343
    dtype: float64




```python
# Compute correlation coefficients of output, inflation, and the nominal interest rate
nk_model.simulated[['y','pi','i']].corr()
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>y</th>
      <th>pi</th>
      <th>i</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>y</th>
      <td>1.000000</td>
      <td>0.697158</td>
      <td>0.756250</td>
    </tr>
    <tr>
      <th>pi</th>
      <td>0.697158</td>
      <td>1.000000</td>
      <td>0.993517</td>
    </tr>
    <tr>
      <th>i</th>
      <td>0.756250</td>
      <td>0.993517</td>
      <td>1.000000</td>
    </tr>
  </tbody>
</table>
</div>


