<jupyter_start><jupyter_text>Getting StartedThis section gives an overview over parsing, simulating and filtering models using **pydsge**. It also explains how to load and process data from an estimation.<jupyter_code># only necessary if you run this in a jupyter notebook
%matplotlib inline

import matplotlib.pyplot as plt
import numpy as np<jupyter_output><empty_output><jupyter_text>Parsing & simulatingLet us first import the base class and the example model:<jupyter_code>from pydsge import * # imports eg. the DSGE class and the examples
yaml_file, data_file = example<jupyter_output><empty_output><jupyter_text>The `example` here is nothing than a tuple containing the paths to two files. The first file to the example model file (as a string):<jupyter_code>print(yaml_file)<jupyter_output>/home/gboehl/repos/pydsge/pydsge/examples/dfi.yaml<jupyter_text>You can use your text editor of choice (which I hope is not Notepad) to open this file and have a look. It contains useful annotations and comments. Parsing ``*.yaml``-filesSo lets parse this thing:<jupyter_code>mod = DSGE.read(yaml_file)<jupyter_output><empty_output><jupyter_text>Of course, if you would want to parse your own ``*.yaml`` model, you could easily do that by defining setting ``yaml_file = "/full/path/to/your/model.yaml"`` instead. But, lets for now assume you are working with ``dfi.yaml``.The ``mod`` object is now an instance of the `DSGE` class. Lets load the calibrated parameters from the file and instantize the transition function:<jupyter_code>par = mod.set_par('calib')<jupyter_output><empty_output><jupyter_text>Have a look at the functions `set_par` and `get_par` in the [modules documentation](https://pydsge.readthedocs.io/en/latest/modules.htmlpydsge.parser.DSGE.get_par) if you want to experiment with the model. `get_par` gets you the current parameters:<jupyter_code># get a single parameter value, e.g. of "kappa"
mod.get_par('kappa')
# or all parameters, for better readability as a dictionary
mod.get_par(asdict=True)<jupyter_output><empty_output><jupyter_text>The ``DSGE``-instance has many attributes that concentrate information and functions of the model. They are all listed in the [module documentation](https://pydsge.readthedocs.io/en/latest/modules.html). Simulate IRFsLet us use this to simulate a series of impulse responses:<jupyter_code>shock_list = ('e_u', 4.0, 0) # (name, size, period)
X1, (L1, K1), _ = mod.irfs(shock_list, verbose=True)<jupyter_output><empty_output><jupyter_text>Nice. For details see the ``irfs`` function. Lets plot it using the ``pplot`` plot function from the ``grgrlib`` library:<jupyter_code>from grgrlib import pplot
figs, axs, _ = pplot(X1)<jupyter_output><empty_output><jupyter_text>Btw, the ``L1, K1`` arrays contain the series of expected durations to/at the ZLB for each simulated period:<jupyter_code>print(L1)
print(K1)<jupyter_output>[0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
 0. 0. 0. 0. 0. 0.]
[7. 6. 5. 4. 3. 2. 1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
 0. 0. 0. 0. 0. 0.]<jupyter_text>Well, lets see what changes if we increase the financial friction parameter `nub`...<jupyter_code>mod.set_par('nub',.5)
X2 = mod.irfs(shock_list)[0] # this is still the same shock as defined above
v = ['Pi','y','r'] # also, lets just look at Pi, y & r
figs, axs, _ = pplot((X1[v],X2[v])) # the output from mod.irfs can be indexed by the variable names!<jupyter_output><empty_output><jupyter_text>Sample from priorNow lets assume that you have specified priors and wanted to know how flexible your model is in terms of impulse responses. The ``get_par`` function also allows sampling from the prior:<jupyter_code>par0 = mod.get_par('prior', nsamples=100, verbose=True)
print(par0.shape)<jupyter_output>(100, 17)<jupyter_text>`par0` is an array with 100 samples of the 10-dimensional parameter vector. If you allow for ``verbose=True`` (which is the default) the function will also tell you how much of your prior is not implicitely trunkated by indetermined or explosive regions. Lets feed these parameters ``par0`` into our ``irfs()`` and plot it:<jupyter_code>mod.set_par('calib')
shock_list = ('e_u', 3., 1) # (name, size, period)

X1, LK1, _ = mod.irfs(shock_list, par0, verbose=True)
v = mod.vix(('y','Pi','x','r')) # get the indices of these four variables
figs, axs, _ = pplot(X1[...,v], labels=mod.vv[v])<jupyter_output><empty_output><jupyter_text>This gives you an idea on how tight your priors are. As the relatively strong negative demand shock kicks the economy to the ZLB, we observe a heavy economic contraction for some of the draws. Filtering & smoothingThis section treats how to load data, and do Bayesian filtering given a DSGE model. Load dataWe have just seen how to parse the model. Parsing the data is likewise quite easy. It however assumes that you managed to put your data into pandas' ``DataFrame`` format. pandas knows many ways of loading your data file into a ``DataFrame``, see for example [here](https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.read_csv.html) on how to load a common ``*.csv`` file. Luckily I already prepared an example data file that is already well structured:<jupyter_code>yaml_file, data_file = example
print(data_file)<jupyter_output>/home/gboehl/repos/pydsge/pydsge/examples/tsdata.csv<jupyter_text>Again, this is just the path to a file that you can open and explore. I constructed the file such that I can already load the column ``data`` as a ``DateTimeIndex``, which makes things easier:<jupyter_code>import pandas as pd

df = pd.read_csv(data_file, parse_dates=['date'], index_col=['date'])
df.index.freq = 'Q' # let pandas know that this is quartely data

print(df)<jupyter_output>GDP     Infl   FFR
date                              
1995-03-31  0.16884  0.60269  1.45
1995-06-30  0.13066  0.44765  1.51
1995-09-30  0.64258  0.42713  1.45
1995-12-31  0.45008  0.47539  1.43
1996-03-31  0.47113  0.50984  1.34
...             ...      ...   ...
2017-03-31  0.41475  0.49969  0.18
2017-06-30  0.54594  0.25245  0.24
2017-09-30  0.54391  0.51972  0.29
2017-12-31  0.48458  0.57830  0.30
2018-03-31  0.15170  0.48097  0.36

[93 rows x 3 columns]<jupyter_text>Now you should give you an idea of how the data looks like. The frame contains the time series of US output growth, inflation, and the FFR from 1995Q1 to 2018Q1. It is generally a good idea to adjust the series of observables such, that they are within reach of the model. As the model can not simulate values below the ZLB, we should adjust the series of the interest rate accordingly:<jupyter_code># adjust elb
zlb = mod.get_par('elb_level')
rate = df['FFR']
df['FFR'] = np.maximum(rate,zlb)
mod.load_data(df)<jupyter_output><empty_output><jupyter_text>This automatically selects the obsevables you defined in the ``*.yaml`` and puts them in the ``mod.data`` object. Note that it will complain if it can't find these observables or if they are named differently. So, that's all we want from now. Run filterWe now want to use a Bayesian Filter to smooth out the hidden states of the model. As the example data sample contains the Zero-lower bound period and the solution method is able to deal with that, we should use a nonlinear filter such as the *Transposed Ensemble Kalman Filter (TEnKF)*. This filter is a hybrid between the Kalman Filter and the Particle Filter, we hence have to define the number of particles. For small problems as the one here, a smaller number would be sufficient, but since everything goes so fast, let us chose 500:<jupyter_code>mod.set_par('calib') # important to reset the parameters as they are still random from previous sampling
mod.create_filter(N=500, ftype='TEnKF', seed=0)<jupyter_output><empty_output><jupyter_text>The **TEnKF** is the default filter, so specifying ``ftype`` would not even have been necessary. The filter got most of the necessary information (innovation covariance, observation function etc) from the ``*.yaml``. What remains to be specified is the measurement noise. The covariance matrix of the measurement errors are stored as ``mod.filter.R``. Luckily, there is a function that creates a diagonal matrix with its diagonal equal to the fraction `a` of the standard deviation of the respective time series, as it is frequently done:<jupyter_code>mod.filter.R = mod.create_obs_cov(1e-1)<jupyter_output><empty_output><jupyter_text>Here, `a=1e-1`. As one last thing before running the filter, we would like to set the ME of the FFR very low as this can be measured directly (note that we can not set it to zero due to numerical reasons, but we can set it sufficiently close).<jupyter_code># lets get the index of the FFR
ind = mod.observables.index('FFR')
# set ME of the FFR to very small value
mod.filter.R[ind,ind] = 1e-4<jupyter_output><empty_output><jupyter_text>``mod.observables`` contains all the observables. See the [module documentation](https://pydsge.readthedocs.io/en/latest/modules.html) for more useful class variables. But lets start the filter already!<jupyter_code>mod.set_par('calib')
FX = mod.run_filter(verbose=True, smoother=True)<jupyter_output>[run_filter:]   Filtering done in 1.055s.<jupyter_text>``smoother=True`` also directly runs the TEnKF-RTS-Smoother. ``FX`` now contains the states. Lets have a look:<jupyter_code>fig, ax, _ = pplot(FX, mod.data.index, labels=mod.vv)<jupyter_output><empty_output><jupyter_text>We can also have a look at the implied observables. The function ``mod.obs()``is the observation function, implemented to work on particle clouds (such as ``FX``):<jupyter_code>FZ = mod.obs(FX)
fig, ax, _ = pplot((mod.Z, FZ), mod.data.index, labels=mod.observables, styles=('-','--'))<jupyter_output><empty_output><jupyter_text>That looks perfect. Note that these particles/ensemble members/"dots" yet do not *fully* obey the nonlinearity of the transition function but contain approximation errors. To get rid of those we need adjustment smoothing. Tools for Model AnalysisThis section shows how to obtain batches of smoothed shocks that can be used to simulate counterfactuals, historic decompositions, and variance decompositions. Extracting/smoothed shocks**pydsge** uses the Nonliear Path-Adjustment Smoother (NPAS) implemented in my [**econsieve**](https://econsieve.readthedocs.io/en/latest/readme.html) package to extract the series of shocks given a transition function and a particle cloud. Keep in mind that in order to account for filtering uncertainty, you need to draw many series from the same particle cloud. For illustrative purposes, lets start with a sample of 20 (which is very likely to be too small):<jupyter_code>mod.set_par('calib') # just to be sure
epd = mod.extract(nsamples=20)<jupyter_output>100%|██████████| 20/20 [00:58<00:00,  2.93s/ sample(s)]<jupyter_text>The `mod.extract` function returns a dictionary (here named `epd`). It contains the parameters of each sample (here, they are all the same), the NPA-smoothed means and covariances, the implied series of observables, and the series of smoothed shock residuals. The `flags` entry contains potential error messages from NPA-smoothing.<jupyter_code>epd.keys()<jupyter_output><empty_output><jupyter_text>Simulating counterfactualsLets use these residuals to simulate counterfactuals. For that, we will use the function `mod.simulate()`, which takes the dictionary of extracted shocks as inputs, potentially together with a mask that will be applied to these shocks to alter them. Let me first create such mask:<jupyter_code>msk0 = mod.mask
print(msk0)<jupyter_output>e_z  e_r  e_u
date                     
1995-03-31  NaN  NaN  NaN
1995-06-30  NaN  NaN  NaN
1995-09-30  NaN  NaN  NaN
1995-12-31  NaN  NaN  NaN
1996-03-31  NaN  NaN  NaN
...         ...  ...  ...
2016-12-31  NaN  NaN  NaN
2017-03-31  NaN  NaN  NaN
2017-06-30  NaN  NaN  NaN
2017-09-30  NaN  NaN  NaN
2017-12-31  NaN  NaN  NaN

[92 rows x 3 columns]<jupyter_text>As you can see this is a dataframe equivalent to the series of shocks, but filled with eps. The logic is simple: whenever the mask is `NaN`, the original shock value for the respective period will be used. If the mask is other then `Nan`, the mask will be used. Let's give this a try by creating a mask that will set all demand shocks (`e_r`) after 2008Q4 equal zero. This can be interpreted as a counterfactual analysis to study the effects of monetary policy shocks and forward guidance during the ZLB period:<jupyter_code>msk0['e_r']['2008Q4':] = 0 
cfs0 = mod.simulate(epd, mask=None)[0] # simulate using original shocks 
cfs1 = mod.simulate(epd, mask=msk0)[0] # simulate using the masked shocks without MP

v = mod.vix(('y','Pi','r','vr'))
fig, ax, _ = pplot((cfs0[...,v], cfs1[...,v]), mod.data.index, labels=mod.vv[v])<jupyter_output><empty_output><jupyter_text>Interesting. According to this toy model, the course of the ZLB is completely independent of monetary policy, yet the policy shocks at the ZLB are able to stimulate inflation through the forward guidance channel.   Historic decompositionsWe can also use the `epd` dictionary to create historic shock decompositions:<jupyter_code>hd, means = mod.nhd(epd)

v = ['y','Pi','r'] # lets only look at these three here

hd = [h[v] for h in hd] # the `hd` object is again a pandas dataframe, so we can index with variable names
hmin, hmax = sort_nhd(hd) # devide into positive and negative contributions

from grgrlib import figurator # the figurator is another nice tool simplifying plotting
fig, ax = figurator(3,1, figsize=(9,7))

pplot(hmax, mod.data.index, labels=v, alpha=.5, ax=ax)
pplot(hmin, mod.data.index, labels=v, alpha=.5, ax=ax, legend=mod.shocks)
pplot(means[v], mod.data.index, labels=v, ax=ax)

fig[0].tight_layout()
_ = ax[2].legend(loc=3)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.0999056809112675, 'fraction_numerical': 0.03584125371834869, 'mean_word_length': 4.801346801346801, 'pattern_counts': {'":': 0, '<': 98, '<?xml version=': 0, '>': 97, 'https://': 5, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '26', 'dir': 'jupyter-structured-clean-dedup', 'id': '73135', 'license': 'MIT', 'n_tokens_mistral': 4801, 'n_tokens_neox': 4332, 'n_words': 1726, 'path': 'docs/getting_started.ipynb', 'repo_name': 'florabudianto/pydsge', 'size': '13783'}