notebook.community

Edit and run 





In [1]:



    


import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import linearsolve as ls
%matplotlib inline

Class 15: Introduction to the linearsolve package

In general, dynamic stochastic general equilibrium (DSGE) models are time-consuming to work with. The linearsolve Python package approximates, solves, and simulates DSGE models.

Example 1: A one-equation model of TFP

Consider the following AR(1) specification for $\log$ TFP:

\begin{align} \log A_{t+1} & = \rho \log A_t + \epsilon_{t+1} \tag{1} \end{align}

where $\epsilon_{t+1} \sim \mathcal{N}(0,\sigma^2)$. Let's simulate the model with linearsolve. To do this we need to do several things:

  1. Create a Pandas series that stores the names of the parameters of the model.
  2. Define a function that returns the equilibrium conditions of the model solved for zero.
  3. Initialize an instance of the linearsolve.model class
  4. Compute and input the steady state of the model.
  5. Approximate and solve the model.
  6. Compute simulations of the model.




In [2]:



    


# 1. Input model parameters
parameters = pd.Series()
parameters['rhoa'] = .9
parameters['sigma'] = 0.001

print(parameters)



    


















    






rhoa     0.900
sigma    0.001
dtype: float64

















In [3]:



    


# 2. Define a function that evaluates the equilibrium conditions
def equilibrium_equations(variables_forward,variables_current,parameters):
    
    # Parameters 
    p = parameters
    
    # Variables
    fwd = variables_forward
    cur = variables_current
        
    # Exogenous tfp
    tfp_proc = p.rhoa*np.log(cur.a) - np.log(fwd.a)
    
    # Stack equilibrium conditions into a numpy array
    return np.array([
        tfp_proc
        ])



    











In [4]:



    


# 3. Initialize the model
model = ls.model(equations = equilibrium_equations,
                 nstates=1,
                 varNames=['a'],
                 shockNames=['eA'],
                 parameters = parameters)



    











In [5]:



    


# 4. Have linearsolve compute the steady state numerically
guess = [1]
model.compute_ss(guess)

print(model.ss)



    


















    






a    1.0
dtype: float64

















In [6]:



    


# 5. Find the log-linear approximation around the non-stochastic steady state and solve
model.approximate_and_solve()



    











In [7]:



    


# 6 (a) Compute impulse responses
model.impulse(T=41,t0=5,shock=None)
print(model.irs['eA'].head(10))



    


















    






          a    eA
0  0.000000  0.00
1  0.000000  0.00
2  0.000000  0.00
3  0.000000  0.00
4  0.000000  0.00
5  0.010000  0.01
6  0.009000  0.00
7  0.008100  0.00
8  0.007290  0.00
9  0.006561  0.00

















In [8]:



    


# 6 (b) Plot impulse responses
model.irs['eA'][['eA','a']].plot(lw='5',alpha=0.5,grid=True).legend(loc='upper right',ncol=2)



    


















    
Out[8]:








<matplotlib.legend.Legend at 0x10ecc4400>










    
























In [9]:



    


# 6(c) Compute stochastic simulation
model.stoch_sim(seed=192,covMat= [[parameters['sigma']]])
print(model.simulated.head(10))



    


















    






          a        eA
0 -0.018781  0.009049
1  0.034836  0.051739
2 -0.035488 -0.066840
3 -0.049489 -0.017550
4 -0.044485  0.000055
5 -0.062945 -0.022909
6 -0.030305  0.026346
7 -0.000027  0.027247
8  0.063010  0.063034
9  0.021703 -0.035006

















In [10]:



    


# 6(d) Plot stochastic simulation
model.simulated[['eA','a']].plot(lw='5',alpha=0.5,grid=True).legend(loc='upper right',ncol=2)



    


















    
Out[10]:








<matplotlib.legend.Legend at 0x10efd2630>

Example 2: A three-equation model business cycle model

Now consider the following system of equations:

\begin{align} Y_t & = A_t K_t^{\alpha} \tag{2}\\ K_{t+1} & = sY_t + (1-\delta) K_t \tag{3}\\ \log A_{t+1} & = \rho \log A_t + \epsilon_{t+1} \tag{4} \end{align}

where $\epsilon_{t+1} \sim \mathcal{N}(0,\sigma^2)$. Let's simulate the model with linearsolve. Before proceding, let's also go ahead and rewrite the model with all variables moved to the lefthand side of the equations:

\begin{align} 0 & = A_t K_t^{\alpha} - Y_t \tag{5}\\ 0 & = sY_t + (1-\delta) K_t - K_{t+1} \tag{6}\\ 0 & = \rho \log A_t + \epsilon_{t+1} - \log A_{t+1} \tag{7} \end{align}

Capital and TFP are called state variables because they're $t+1$ values are predetermined. Output is called a costate variable. Note that the model as 3 endogenous variables with 2 state variables.





In [11]:



    


# 1. Input model parameters
parameters = pd.Series()
parameters['rhoa'] = .9
parameters['sigma'] = 0.01
parameters['alpha'] = 0.35
parameters['delta'] = 0.025
parameters['s'] = 0.15

print(parameters)



    


















    






rhoa     0.900
sigma    0.010
alpha    0.350
delta    0.025
s        0.150
dtype: float64

















In [12]:



    


# 2. Define a function that evaluates the equilibrium conditions
def equilibrium_equations(variables_forward,variables_current,parameters):
    
    # Parameters 
    p = parameters
    
    # Variables
    fwd = variables_forward
    cur = variables_current
        
    # Production function
    prod_fn = cur.a*cur.k**p.alpha - cur.y
    
    # Capital evolution
    capital_evolution = p.s*cur.a*cur.k**p.alpha + (1 - p.delta)*cur.k - fwd.k
    
    # Exogenous tfp
    tfp_proc = p.rhoa*np.log(cur.a) - np.log(fwd.a)
    
    
    # Stack equilibrium conditions into a numpy array
    return np.array([
        prod_fn,
        capital_evolution,
        tfp_proc
        ])



    











In [13]:



    


# 3. Initialize the model
model = ls.model(equations = equilibrium_equations,
                 nstates=2,
                 varNames=['a','k','y'],                   # Any order as long as the state variables are named first
                 shockNames=['eA','eK'],                 # Name a shock for each state variable *even if there is no corresponding shock in the model*
                 parameters = parameters)



    











In [14]:



    


# 4. Have linearsolve compute the steady state numerically
guess = [1,4,1]
model.compute_ss(guess)



    











In [15]:



    


# 5. Find the log-linear approximation around the non-stochastic steady state and solve
model.approximate_and_solve()

The previous step constructs a log-linear approximation of the model and then solves for the endogenous variables as functions of the state variables and exogenous shocks only:

\begin{align} \left[ \hat{y}_t\right] & = F \left[\begin{array}{c} \hat{a}_t \\ \hat{k}_t \end{array}\right]\\ \left[\begin{array}{c} \hat{a}_{t+1} \\ \hat{k}_{t+1} \end{array}\right] & = P \left[\begin{array}{c} \hat{a}_t \\ \hat{k}_t \end{array}\right] + \left[\begin{array}{c} \epsilon_{t+1}^a \\ \epsilon^k_{t+1} \end{array}\right]. \end{align}

where $F$ and $P$ are coefficient matrices computed by the program.





In [16]:



    


# Print the coeficient matrix P
print(model.p)



    


















    






[[ 0.9      0.     ]
 [ 0.025    0.98375]]

















In [17]:



    


# Print the coeficient matrix F
print(model.f)



    


















    






[[ 1.    0.35]]

















In [18]:



    


# 6 (a) Compute impulse responses and print the computed impulse responses
model.impulse(T=41,t0=5,shock=None)
print(model.irs['eA'].head(10))



    


















    






          a    eA         k         y
0  0.000000  0.00  0.000000  0.000000
1  0.000000  0.00  0.000000  0.000000
2  0.000000  0.00  0.000000  0.000000
3  0.000000  0.00  0.000000  0.000000
4  0.000000  0.00  0.000000  0.000000
5  0.010000  0.01  0.000000  0.010000
6  0.009000  0.00  0.000250  0.009088
7  0.008100  0.00  0.000471  0.008265
8  0.007290  0.00  0.000666  0.007523
9  0.006561  0.00  0.000837  0.006854

















In [19]:



    


# 6(b) Plot the computed impulse responses to a TFP shock
fig = plt.figure(figsize=(12,4))

ax1 = fig.add_subplot(1,2,1)
model.irs['eA'][['a','y','k']].plot(lw='5',alpha=0.5,grid=True,ax = ax1).legend(loc='upper right',ncol=2)

ax2 = fig.add_subplot(1,2,2)
model.irs['eA'][['eA','a']].plot(lw='5',alpha=0.5,grid=True,ax = ax2).legend(loc='upper right',ncol=2)



    


















    
Out[19]:








<matplotlib.legend.Legend at 0x10f33e390>










    
























In [20]:



    


# 6(c) Compute stochastic simulation and print the simulated values
model.stoch_sim(seed=192,covMat= [[parameters['sigma'],0],[0,0]])
print(model.simulated.head(10))



    


















    






          a        eA   eK         k         y
0 -0.045037  0.014046  0.0 -0.227059 -0.124508
1 -0.054616 -0.014082  0.0 -0.224495 -0.133189
2 -0.090089 -0.040935  0.0 -0.222213 -0.167863
3 -0.139656 -0.058576  0.0 -0.220854 -0.216955
4 -0.186526 -0.060835  0.0 -0.220756 -0.263791
5 -0.216183 -0.048309  0.0 -0.221832 -0.293824
6 -0.208871 -0.014307  0.0 -0.223632 -0.287142
7 -0.090173  0.097811  0.0 -0.225220 -0.169000
8 -0.067672  0.013483  0.0 -0.223814 -0.146007
9 -0.181917 -0.121012  0.0 -0.221869 -0.259571

















In [21]:



    


# 6(d) Plot the computed stochastic simulation
fig = plt.figure(figsize=(12,4))

ax1 = fig.add_subplot(1,2,1)
model.simulated[['a','y','k']].plot(lw='5',alpha=0.5,grid=True,ax = ax1).legend(loc='upper right',ncol=3)

ax2 = fig.add_subplot(1,2,2)
model.simulated[['eA','a']].plot(lw='5',alpha=0.5,grid=True,ax = ax2).legend(loc='upper right',ncol=2)



    


















    
Out[21]:








<matplotlib.legend.Legend at 0x10f834fd0>

Content source: letsgoexploring/teaching

Similar notebooks:

notebook.community | gallery | about