# PyDSGE Example - Small New Keynesian Model
## Author: Gustavo Amarante
This jupyter notebook is an example on how to implement the basic new keynesian model. The details of the model and 
parameters interpretations can be found in [Galí (2008)](https://www.amazon.com/Monetary-Policy-Inflation-Business-Cycle/dp/0691133166) 
and a full step-by-step derivation of the model equations is in 
[Bergholt (2012)](https://bergholt.weebly.com/uploads/1/1/8/4/11843961/the_basic_new_keynesian_model_-_drago_bergholt.pdf).

The linearized version of the model is given by the following equations:

$
\begin{align*}
\tilde{y}_{t} & =E_{t}\left(\tilde{y}_{t+1}\right)-\frac{1}{\sigma}\left[\hat{i}_{t}-E_{t}\left(\pi_{t+1}\right)\right]+\psi_{ya}^{n}\left(\rho_{a}-1\right)a_{t}\\
\pi_{t} & =\beta E_{t}\left(\pi_{t+1}\right)+\kappa\tilde{y}_{t}+\sigma_{\pi}\varepsilon_{t}^{\pi}\\
\hat{i}_{t} & =\phi_{\pi}\pi_{t}+\phi_{y}\tilde{y}_{t}+v_{t}\\
a_{t} & =\rho_{a}a_{t-1}+\sigma_{a}\varepsilon_{t}^{a}\\
v_{t} & =\rho_{v}v_{t-1}+\sigma_{v}\varepsilon_{t}^{v}
\end{align*}
$

where the following parameters are given by:

$
\psi_{ya}^{n}=\frac{1+\varphi}{\sigma\left(1-\alpha\right)+\varphi+\alpha}\qquad\kappa=\frac{\left(1-\theta\right)\left(1-\theta\beta\right)\left[\sigma\left(1-\alpha\right)+\varphi+\alpha\right]}{\theta\left(1-\alpha+\alpha\epsilon\right)}
$

We also need to write the expectations of time $t+1$ variables as time-$t$ variables. The relationship between variable 
expectations and the effective value is just the expectational error.

$
\begin{align*}
\tilde{y}_{t+1} & =ex_{t}^{y}+\eta_{t+1}^{y}\\
\pi_{t+1} & =ex_{t}^{\pi}+\eta_{t+1}^{\pi}
\end{align*}
$

where $ex_{t}^{y}=E_{t}\left(\tilde{y}_{t+1}\right)$.

So we have 7 endogenous variables, 3 exogenous shocks $\varepsilon$ and 2 expectational errors $\eta$ that compose the 
state equations of the model. We also need to specify the set of observation equations, which we will assume that only 
the outputgap, inflation and interest rate are observable variables. 

$
\begin{align*}
output\_gap_{t} & =\tilde{y}_{t}\\
inflation_{t} & =\pi_{t}\\
interest\_rate_{t} & =\left(\frac{1}{\beta}-1\right)+\hat{i}_{t}\\
\end{align*}
$

These are the sections of the notebook:
* Especification
* Calibration
* Simulation
* Estimation
* Analysis

---
# Imports
Your imports should include `symbols` and `Matrix` from the [`sympy` library](https://www.sympy.org/en/index.html), as 
these are required for variables and model declaration.

(Obs: the need for libraries should be eliminated in future versions) 
 


```python
import numpy as np
from pydsge import DSGE
import matplotlib.pyplot as plt
from sympy import symbols, Matrix
```

---
# Model Especification
## Variable Declaration
Declare all of the endogenous variables at time-$t$ as symbols and create a `Matrix` object with them.


```python
y, pi, i, a, v, exp_y, exp_pi = symbols('y, pi, i, a, v, exp_y, exp_pi')
endog = Matrix([y, pi, i, a, v, exp_y, exp_pi])
```

Declare all of the endogenous variables at time-$t-1$ as symbols and create a `Matrix` object with them. The convention 
in this notebook is to  add an `l` (for "lagged") at the end of the variable name.


```python
yl, pil, il, al, vl, exp_yl, exp_pil = symbols('yl, pil, il, al, vl, exp_yl, exp_pil')
endogl = Matrix([yl, pil, il, al, vl, exp_yl, exp_pil])
```

Declare all of the exogenous shocks $\varepsilon$ as symbols and create a `Matrix` object with them.


```python
eps_a, eps_v, eps_pi = symbols('eps_a, eps_v, eps_pi')
exog = Matrix([eps_a, eps_v, eps_pi])
```

Declare all of the expectational errors $\eta$ as symbols and create a `Matrix` object with them.


```python
eta_y, eta_pi = symbols('eta_y, eta_pi')
expec = Matrix([eta_y, eta_pi])

```

Declare all of the parameter as symbols and create a `Matrix` object with them. 
Summary parameters (functions of other parameters) do not need to be inside the matrix object, they only need to be 
defined by their formula.


```python
sigma, varphi, alpha, beta, theta, phi_pi, phi_y, rho_a, sigma_a, rho_v, sigma_v, sigma_pi = \
    symbols('sigma, varphi, alpha, beta, theta, phi_pi, phi_y, rho_a, sigma_a, rho_v, sigma_v, sigma_pi')

psi_nya = (1 + varphi) / (sigma*(1-alpha) + varphi + alpha)
kappa = (1 - theta)*(1 - theta * beta)*(sigma*(1-alpha) + varphi + alpha)
```

## State Equations 
Given the set of linearized equilibrium conditions described in the beggining of thie notebook, put all of the terms 
to one side of the equation so that the system of equations is equal to zero. Since all of our variable and parameter 
names are symbols, the equations will be correctly interpreted as symbolic expressions. These expressions should also 
go into a `Matrix` object. This set of equations should include only the state equations of the system. 


```python
eq1 = y - exp_y + (1/sigma)*(i - exp_pi) - psi_nya * (rho_a - 1) * a
eq2 = pi - beta * exp_pi - kappa * y - sigma_pi * eps_pi
eq3 = i - phi_pi * pi - phi_y * y - v
eq4 = a - rho_a * al - sigma_a * eps_a
eq5 = v - rho_v * vl - sigma_v * eps_v
eq6 = y - exp_yl - eta_y
eq7 = pi - exp_pil - eta_pi

equations = Matrix([eq1, eq2, eq3, eq4, eq5, eq6, eq7])
```

## Observation Equations
Using the observation equations described in the beggining of thie notebook, the expression for the observable variables
(linear function of the state variables) should be on the right side of the equation. These expressions should also go 
into a `Matrix` object. This set of equations  should include only the observation equations of the system. 


```python
obs01 = y
obs02 = pi
obs03 = 1/beta - 1 + i

obs_equations = Matrix([obs01, obs02, obs03])
```

---
# Model Calibration
Now that we have decribed the full set of equilibrium conditions of the models, we can calibrate the parameters. We will
save the simulated data to illuestrate the estimation processo later on. To calibrate a model, all you need is a 
dictionary where the keys are the parameter symbols (*not* the string with their names) and the values are their 
respective calibration value. The values below will be our "true" parameters. Later in this notebook we will estimate
them.


```python
calib_dict = {sigma: 1.3,
              varphi: 1,
              alpha: 0.4,
              beta: 0.997805,
              theta: 0.75,
              phi_pi: 1.5,
              phi_y: 0.2,
              rho_a: 0.9,
              sigma_a: 1.1,
              rho_v: 0.5,
              sigma_v: 0.3,
              sigma_pi: 0.8}
```

to build the DSGE object we just have to pass the `Matrix` objects that we created and the calibration dictionary.


```python
dsge_simul = DSGE(endog=endog, 
                  endogl=endogl, 
                  exog=exog, 
                  expec=expec, 
                  state_equations=equations,
                  obs_equations=obs_equations,
                  calib_dict=calib_dict)
```

Since we passed a calibration dictionary, the model already knows that it only has to substitute the calibration values
in their respective symbolic representaions and solves the model. Just by writing the statement above, the `dsge_simul`
object already computes everything.

To check if the calbration values yield a solution or not, the model has the `.eu` attibute, mirrors the definition from
the original *Sims (2002)* `gensys.m` function. It represents the existance and uniqueness of the solaution of the model. 
If the first entry is equal to 1 it means the model has a solution and if the second entry is equal to 1 it means that
the solution is unique.


```python
print(dsge_simul.eu)
```

    [1, 1]


---
# Simulation

There is a method to simulate data from the calibrated model. Here we are using a random seed for replication purposes.
We will generate the data and plot the observed variables.


```python
df_obs, df_states = dsge_simul.simulate(n_obs=250, random_seed=1)

df_obs.plot(figsize=(15, 9))
```

---
# Estimation
The estimation DSGEs usually comes with the problem of identificiationm of some of their parameters. To handle this 
problem, the parameters that are not identified usually are calibrated prior to estimation. In the cell below, we first 
create a dictionary of the parameters we want to estimate and then create a `Matrix` object of the parameters that are 
going to be estimated.


```python
calib_param = {varphi: 1, alpha: 0.4, beta: 0.997805}
estimate_param = Matrix([sigma, theta, phi_pi, phi_y, rho_a, sigma_a, rho_v, sigma_v, sigma_pi])
```

Since we use bayesian estimation with MCMC we need to define the priors of the parameters that are going to be 
estimated. The available distributions for priors are Beta, Gamma, Inverse Gamma, Uniform and Normal. For each prior you 
have to choose a mean and a standard deviation. The table below shows how the mean and standard deviations are related 
to the distribution parameters. There is also a `label` entry for each parameter, which will be used for legends and 
titles on the output charts, and notice that it accepts $\LaTeX$ entries.

|Distribution|Parameter $a$|Parameter $b$|Support           |Mean           |Variance                                           |
|------------|-------------|-------------|------------------|---------------|---------------------------------------------------|
|Beta        |$a>0$ (shape)|$b>0$ (shape)|$\left[0,1\right]$|$\frac{a}{a+b}$|$\frac{ab}{\left(a+b\right)^{2}\left(a+b+1\right)}$|
|Gamma	     |$a>0$ (shape)|$b>0$ (scale)|$\mathbb{R}^{+}$  |$ab$           |$ab^{2}$                                           |
|Inverse Gamma|$a>0$ (shape)|$b>0$ (scale)|$\mathbb{R}^{+}$ |$\frac{\beta}{\alpha-1}$|$\frac{\beta^{2}}{\left(\alpha-1\right)^{2}\left(\alpha-2\right)}$|
|Uniform     |$a\in\mathbb{R}$ (lower bound)|b>a (upper bound)|$\left[a,b\right]$|$\frac{a+b}{2}$|$\frac{\left(b-a\right)^{2}}{12}$|
|Normal      |$a\in\mathbb{R}$ (mean)|$b>0$ (std)|$\mathbb{R}$|$a$|$b^{2}$|


```python
prior_dict = {sigma:    {'dist': 'normal',   'mean':  1.30, 'std': 0.20, 'label': '$\\sigma$'},
              theta:    {'dist': 'beta',     'mean':  0.60, 'std': 0.20, 'label': '$\\theta$'},
              phi_pi:   {'dist': 'normal',   'mean':  1.50, 'std': 0.35, 'label': '$\\phi_{\\pi}$'},
              phi_y:    {'dist': 'gamma',    'mean':  0.25, 'std': 0.10, 'label': '$\\phi_{y}$'},
              rho_a:    {'dist': 'beta',     'mean':  0.50, 'std': 0.25, 'label': '$\\rho_a$'},
              sigma_a:  {'dist': 'invgamma', 'mean':  0.50, 'std': 0.25, 'label': '$\\sigma_a$'},
              rho_v:    {'dist': 'beta',     'mean':  0.50, 'std': 0.25, 'label': '$\\rho_v$'},
              sigma_v:  {'dist': 'invgamma', 'mean':  0.50, 'std': 0.25, 'label': '$\\sigma_v$'},
              sigma_pi: {'dist': 'invgamma', 'mean':  0.50, 'std': 0.25, 'label': '$\\sigma_{\\pi}$'}}

dsge = DSGE(endog=endog, 
            endogl=endogl, 
            exog=exog, 
            expec=expec, 
            state_equations=equations,
            estimate_params=estimate_param,
            calib_dict=calib_param,
            obs_equations=obs_equations,
            prior_dict=prior_dict,
            obs_data=df_obs,
            verbose=True)
```

After the model is declared, you can to pass three arguments to the `estimate` method. 

The argument `file_path` is the address for where you would like to save the HDF file with the generated MCMC chains. 
This filetype was chosen because of its fast read/write speed (on the other hand, the file is very large). Every 100 
iterations, the chains will be saved to this file, which ensure that the MCMC chains are not lost on the event that the 
execution is interrupted. If you pass the adress of a file that already exists, it will read the file and continue the 
chain from where it stoped.

The `nsim` (default: 1000) is the size of the chain and `ck` (deafult: 0.2) is the scale parameter for the
Accpetance-Rejection part of the metropolis-hasting algorithm, which should be set in order to match your desired 
acceptance rate.


```python
dsge.estimate(file_path='example_snkm.h5', nsim=2000, ck=0.3)
```

    Metropolis-Hastings:   3%|▎         | 63/1999 [00:36<17:41,  1.82it/s]/Users/gustavoamarante/anaconda3/envs/py36/lib/python3.6/site-packages/ipykernel_launcher.py:1: PerformanceWarning: 
    your performance may suffer as PyTables will pickle object types that it cannot
    map directly to c-types [inferred_type->floating,key->block0_values] [items->Index(['sigma', 'theta', 'phi_pi', 'phi_y', 'rho_a', 'sigma_a', 'rho_v',
           'sigma_v', 'sigma_pi'],
          dtype='object')]
    
      """Entry point for launching an IPython kernel.
    Metropolis-Hastings: 100%|██████████| 1999/1999 [18:13<00:00,  1.83it/s]


    Acceptance rate: 28.349999999999998 percent


The `eval_chain` method generates two charts. The first one contains the subplots of the MCMC chains and the second one
is the prior (countinous line) and posterior (histogram) for each parameter. The `burnin` parameter defines how much of
the beggining of the chain left out of the analysis. If the input is a `float` between 0 and 1, that percentage of the
beggining of the chain will be dropped, and if it is an `int`, that many observations are dropped.


```python
dsge.eval_chains(burnin=0, show_charts=True)
```

After the estimation, the DSGE model obsject earns the `posterior_table` attribute, which contains a pandas DataFrame
with the tradional measures reported on academic papers.


```python
dsge.posterior_table
```




<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>prior dist</th>
      <th>prior mean</th>
      <th>prior std</th>
      <th>posterior mode</th>
      <th>posterior mean</th>
      <th>posterior 5%</th>
      <th>posterior 95%</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>sigma</th>
      <td>normal</td>
      <td>1.3</td>
      <td>0.2</td>
      <td>1.225852</td>
      <td>1.235096</td>
      <td>1.226802</td>
      <td>1.242856</td>
    </tr>
    <tr>
      <th>theta</th>
      <td>beta</td>
      <td>0.6</td>
      <td>0.2</td>
      <td>0.673278</td>
      <td>0.678953</td>
      <td>0.674504</td>
      <td>0.682882</td>
    </tr>
    <tr>
      <th>phi_pi</th>
      <td>normal</td>
      <td>1.5</td>
      <td>0.35</td>
      <td>1.496686</td>
      <td>1.500633</td>
      <td>1.489683</td>
      <td>1.512747</td>
    </tr>
    <tr>
      <th>phi_y</th>
      <td>gamma</td>
      <td>0.25</td>
      <td>0.1</td>
      <td>0.177505</td>
      <td>0.175000</td>
      <td>0.171605</td>
      <td>0.178876</td>
    </tr>
    <tr>
      <th>rho_a</th>
      <td>beta</td>
      <td>0.5</td>
      <td>0.25</td>
      <td>0.744787</td>
      <td>0.747697</td>
      <td>0.739498</td>
      <td>0.754306</td>
    </tr>
    <tr>
      <th>sigma_a</th>
      <td>invgamma</td>
      <td>0.5</td>
      <td>0.25</td>
      <td>0.529376</td>
      <td>0.530659</td>
      <td>0.524724</td>
      <td>0.536755</td>
    </tr>
    <tr>
      <th>rho_v</th>
      <td>beta</td>
      <td>0.5</td>
      <td>0.25</td>
      <td>0.377100</td>
      <td>0.377956</td>
      <td>0.371537</td>
      <td>0.384737</td>
    </tr>
    <tr>
      <th>sigma_v</th>
      <td>invgamma</td>
      <td>0.5</td>
      <td>0.25</td>
      <td>0.529523</td>
      <td>0.533092</td>
      <td>0.527101</td>
      <td>0.539630</td>
    </tr>
    <tr>
      <th>sigma_pi</th>
      <td>invgamma</td>
      <td>0.5</td>
      <td>0.25</td>
      <td>0.675716</td>
      <td>0.668371</td>
      <td>0.663087</td>
      <td>0.673798</td>
    </tr>
  </tbody>
</table>
</div>



ANALIZAR A TABELA DAS POSTERIORIS
REVISAR COMO AS CADEIAS SÃO FORMADAS

618832e07386912ddb307e9fd8947e5e612e91b8
ipynb
Jupyter Notebook
Example/example_snkm.ipynb
gusamarante/pydsge
f018b38a95bf65aaf2a54b64b3140c56d863c61e
MIT
2020-08-03T23:57:07.000Z
2022-03-09T21:28:06.000Z
Example/example_snkm.ipynb
gusamarante/pydsge
f018b38a95bf65aaf2a54b64b3140c56d863c61e
MIT
Example/example_snkm.ipynb
gusamarante/pydsge
f018b38a95bf65aaf2a54b64b3140c56d863c61e
MIT
2020-12-07T04:27:44.000Z
2022-01-31T14:39:38.000Z
Qwen/Qwen-72B
1. YES
2. YES