Estimation Tutorial#
In this section, we dive into the topic of model estimation using pydsge.
Let us, just for the sake of this tutorial, set up a temporary directory structure:
[1]:
# Just for the tutorial: Setting up example structure
import tempfile
import os
import shutil # For clean-up of temporary directory
from pathlib import Path # For Windows/Unix compatibility
# Temporary output folder
output_path = Path(tempfile.gettempdir(), 'output')
if not os.path.isdir(output_path):
os.makedirs(output_path)
Parsing and loading the model#
Let us first load the relevant packages. Besides the DSGE class we already know from getting started, we also want to import the emcee package. This will allow us to later specify the desired updating algorithms for sampling from the posterior distribution - we explain this in more detail below.
[2]:
import pandas as pd
import numpy as np
import emcee # For specifying updating moves
from pydsge import DSGE, example
In this tutorial, we continue to use the example provided in pydsge. Like before, we specify the file paths of the model and the data. Please feel free to check-out both files, but from the previous tutorial you might remember that we’re dealing with a five equations New Keynesian model and US quarterly data from 1995 to 2018.
[3]:
yaml_file, data_file = example
We again parse the model and load the data. Importantly we also specify a location where the output is stored.
For this tutoral I assign a tempfile as the output path, but this is certainly not what you want. So better change that path. Note also that we can give the model estimation a name and also write a short description, which is very useful when working with several models or estimations.
[4]:
# Parse the model
mod = DSGE.read(yaml_file)
# Give it a name
mod.name = 'rank_tutorial'
mod.description = 'RANK, estimation tutorial'
# Storage location for output. CHANGE THIS to some local path on your PC!
mod.path = output_path
# Load data
df = pd.read_csv(data_file, parse_dates=['date'], index_col=['date'])
df.index.freq = 'Q' # let pandas know that this is quartely data
Remember that since the Great Recession, the Federal Funds Rate has been below the ZLB. That is why, like in getting started, we adjust the observed interest rate, so that the data is “within reach” of our model.
[5]:
# adjust elb
zlb = mod.get_par('elb_level')
rate = df['FFR']
df['FFR'] = np.maximum(rate,zlb)
mod.load_data(df, start='1998Q1')
[5]:
| GDP | Infl | FFR | |
|---|---|---|---|
| date | |||
| 1998-03-31 | 0.77834 | 0.14386 | 1.38 |
| 1998-06-30 | 0.69635 | 0.22873 | 1.38 |
| 1998-09-30 | 1.03077 | 0.36109 | 1.38 |
| 1998-12-31 | 1.37921 | 0.26145 | 1.22 |
| 1999-03-31 | 0.54307 | 0.37393 | 1.18 |
| ... | ... | ... | ... |
| 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 |
81 rows × 3 columns
Preparing the estimation#
After importing the packages and loading the data, we still need to tell pydsge how to carry out the estimation of our model. The “prep_estim” method can be used to accomplish this. It can be called without any arguments and sets-up a non-linear model by default. However, not all defaults are always a good good choice, and to showcase some of this functionality, we decide to specify several arguments here.
To perform the estimation, pydsge uses a Transposed-Ensemble Kalman Filter (TEnKF). For general information on its implementation, see the EconSieve documentation , and for more details on running the filter in pydsge check-out the getting started tutorial. Again, the default filter is non-linear, but we can opt for a linear one by setting the argument linear to
True. To choose a custom number of ensemble members for the TEnKF, set N to a particular number (default is 300, for e.g. a medium scale model 400-500 is a good choice). We can also set a specific random seed with the argument seed (the default seed is 0). To get additional information on the estimation process, we can set verbose to True. Conveniently, this information includes an overview of the parameters’ distribution, their means and standard deviations. Finally, if
we already specified the covariance matrix of the measurement errors or want to reuse a previous result, we can load it into the prep_estim method by setting Load.R to True.
If you run into problems you can turn parallelization off by setting debug=True.
[6]:
mod.prep_estim(N=350, seed=0, verbose=True, use_prior_transform=True)
[estimation:] Model operational. 12 states, 3 observables, 3 shocks, 81 data points.
Adding parameters to the prior distribution...
- theta as beta (0.5, 0.1). Init @ 0.7813, with bounds (0.2, 0.95)
- sigma as normal (1.5, 0.375). Init @ 1.2312, with bounds (0.25, 3)
- phi_pi as normal (1.5, 0.25). Init @ 1.7985, with bounds (1.0, 3)
- phi_y as normal (0.125, 0.05). Init @ 0.0893, with bounds (0.001, 0.5)
- rho_u as beta (0.5, 0.2). Init @ 0.7, with bounds (0.01, 0.9999)
- rho_r as beta (0.5, 0.2). Init @ 0.7, with bounds (0.01, 0.9999)
- rho_z as beta (0.5, 0.2). Init @ 0.7, with bounds (0.01, 0.9999)
- rho as beta (0.75, 0.1). Init @ 0.8, with bounds (0.5, 0.975)
- sig_u as inv_gamma_dynare (0.1, 2). Init @ 0.5, with bounds (0.025, 5)
- sig_r as inv_gamma_dynare (0.1, 2). Init @ 0.5, with bounds (0.01, 3)
- sig_z as inv_gamma_dynare (0.1, 2). Init @ 0.5, with bounds (0.01, 3)
[estimation:] 11 priors detected. Adding parameters to the prior distribution.
You probably want to set use_prior_transform to True, as above. This activates the destinction between prior and sampling space for the ADEMC sampler, which reduces the number of neccessary MCMC draws a lot.
As in the filtering tutorial, we set the covariance of measurement errors to correspond to the variances of the data. Additionally, we adjust the measurement errors of the Federal Funds rate since it is perfectly observable.
[7]:
mod.filter.R = mod.create_obs_cov(1e-1)
ind = mod.observables.index('FFR')
mod.filter.R[ind,ind] /= 1e1
Running the estimation#
Lets turn to the actual estimation. For a variety of pretty good reasons, pdygse uses Ensemble Markov Chain Monte Carlo (Ensemble-MCMC) integration to sample from the posterior distribution. For further information on Ensemble-MCMC, please refer to the emcee website and the additional resources provided there.
We first require an initial ensemble, which is provided by tmcmc. tmcmc is a very sophisticated function with many options, but right now, all we are interested in is to obtain a sample that represents the prior distribution:
[8]:
p0 = mod.prior_sampler(40, verbose=True) # rule of thumb: number_of_parameters times 4
100%|██████████| 40/40 [00:06<00:00, 5.83it/s]
(prior_sample:) Sampling done. Check fails for 4.76% of the prior.
The parameter draws are saved in the object p0 as a numpy array in order to later pass them to our main sampling process.
[9]:
mod.save()
[save_meta:] Metadata saved as '/tmp/output/rank_tutorial_meta'
mod.save() saved the meta data of our model in the directory which we specified earlier in mod.path. This information is stored as an .npz file so that it is avialable even in the event of a crash and can be loaded anytime using numpy.load().
For posterior sampling using mcmc we have the option to set different “moves”, i.e. coordinate updating algorithms for the walkers. As a wrapper for a lot of emcee functionality, mcmc can work with many different “moves” - for a list and implementation details please consult the emcee documentation.
The default is however to use the DIME MCMC sampler from this paper, which is very efficient for burn-in and robust to multimodal distributions. The sampler is provided by the emcwrap package. This is also the default sampler for pydsge, so the below is actually not necessary and the moves keyword can just be omited.
[10]:
from emcwrap import DIMEMove
moves = DIMEMove(aimh_prob=0.1)
We now use the initial states derived above to conduct our full Bayesian estimation using mcmc. Note that, instead of using the specified initial ensemble, mcmc can identify previous runs or estimations, or the initial values of the “prior” section in the *.yaml can be used.
The default number of sampling steps is 3000, which is parallelized by default. With tune we can determine the size of the Markov Chain we wish to retain to represent the posterior, i.e. after burn-in. This is not to be confused this with the updating frequency, which only affects the number of summary statements pydsgereports during the estimation.
With the option lprob_seed the user can choose how to set the random seed of the likelihood evaluation - here we use the seed specified in prep_estim.
[11]:
mod.mcmc(p0,
moves=moves,
nsteps=400,
tune=100,
update_freq=50,
) # this may take some time. Better run on a machine with MANY cores...
mod.save() # be sure to save the internal state!
[mcmc:] HDF backend at /tmp/output/rank_tutorial_sampler.h5 already exists. Deleting...
[ll/MAF:-124.691(3e+01)/10% | -2e+01]: 13%|█▎ | 51/400 [03:12<20:09, 3.47s/sample(s)]
(mcmc:) Summary from last 50 of 50 iterations (RANK, estimation tutorial):
distribution pst_mean sd/df mean sd mode hpd_5 hpd_95 \
theta beta 0.500 0.100 0.606 0.067 0.625 0.493 0.697
sigma normal 1.500 0.375 1.108 0.537 1.446 0.246 1.864
phi_pi normal 1.500 0.250 1.821 0.640 1.521 0.993 2.891
phi_y normal 0.125 0.050 0.135 0.156 0.134 -0.129 0.380
rho_u beta 0.500 0.200 0.725 0.072 0.610 0.610 0.835
rho_r beta 0.500 0.200 0.636 0.100 0.663 0.517 0.803
rho_z beta 0.500 0.200 0.420 0.292 0.576 0.017 0.847
rho beta 0.750 0.100 0.724 0.137 0.686 0.520 0.926
sig_u inv_gamma_dynare 0.100 2.000 0.757 0.274 1.067 0.300 1.139
sig_r inv_gamma_dynare 0.100 2.000 0.666 0.459 1.118 0.086 1.221
sig_z inv_gamma_dynare 0.100 2.000 1.063 0.316 1.053 0.563 1.505
error
theta 0.006
sigma 0.042
phi_pi 0.061
phi_y 0.016
rho_u 0.007
rho_r 0.010
rho_z 0.029
rho 0.014
sig_u 0.015
sig_r 0.039
sig_z 0.030
Mean acceptance fraction: 0.249
Autocorrelation times are between 2.64 and 3.56.
[ll/MAF:-59.472(2e+01)/ 5% | -5e-01]: 25%|██▌ | 101/400 [06:17<20:01, 4.02s/sample(s)]
(mcmc:) Summary from last 50 of 100 iterations (RANK, estimation tutorial):
distribution pst_mean sd/df mean sd mode hpd_5 hpd_95 \
theta beta 0.500 0.100 0.558 0.166 0.584 0.314 0.814
sigma normal 1.500 0.375 1.135 0.464 1.420 0.368 1.833
phi_pi normal 1.500 0.250 2.668 0.929 2.983 1.024 3.758
phi_y normal 0.125 0.050 0.077 0.197 -0.138 -0.212 0.379
rho_u beta 0.500 0.200 0.771 0.078 0.747 0.639 0.873
rho_r beta 0.500 0.200 0.570 0.161 0.556 0.322 0.822
rho_z beta 0.500 0.200 0.663 0.414 0.936 0.016 0.987
rho beta 0.750 0.100 0.706 0.165 0.544 0.432 0.917
sig_u inv_gamma_dynare 0.100 2.000 0.384 0.148 0.515 0.134 0.573
sig_r inv_gamma_dynare 0.100 2.000 0.324 0.221 0.673 0.063 0.595
sig_z inv_gamma_dynare 0.100 2.000 0.888 0.345 0.727 0.378 1.505
error
theta 0.017
sigma 0.036
phi_pi 0.097
phi_y 0.020
rho_u 0.008
rho_r 0.017
rho_z 0.042
rho 0.018
sig_u 0.014
sig_r 0.023
sig_z 0.028
Mean acceptance fraction: 0.208
Autocorrelation times are between 4.78 and 7.75.
[ll/MAF:-18.338(9e+00)/12% | -8e-02]: 38%|███▊ | 151/400 [09:46<17:47, 4.29s/sample(s)]
(mcmc:) Summary from last 50 of 150 iterations (RANK, estimation tutorial):
distribution pst_mean sd/df mean sd mode hpd_5 hpd_95 \
theta beta 0.500 0.100 0.710 0.183 0.947 0.480 0.981
sigma normal 1.500 0.375 1.735 0.579 1.702 0.967 2.833
phi_pi normal 1.500 0.250 2.743 0.895 0.561 1.228 4.030
phi_y normal 0.125 0.050 0.088 0.153 0.531 -0.113 0.360
rho_u beta 0.500 0.200 0.897 0.047 0.851 0.837 0.970
rho_r beta 0.500 0.200 0.575 0.196 0.878 0.274 0.871
rho_z beta 0.500 0.200 0.766 0.374 0.002 0.018 0.999
rho beta 0.750 0.100 0.837 0.077 0.814 0.727 0.934
sig_u inv_gamma_dynare 0.100 2.000 0.252 0.120 0.649 0.095 0.395
sig_r inv_gamma_dynare 0.100 2.000 0.169 0.088 0.107 0.061 0.269
sig_z inv_gamma_dynare 0.100 2.000 0.406 0.374 0.424 0.047 0.842
error
theta 0.020
sigma 0.064
phi_pi 0.109
phi_y 0.019
rho_u 0.004
rho_r 0.023
rho_z 0.046
rho 0.009
sig_u 0.013
sig_r 0.010
sig_z 0.040
Mean acceptance fraction: 0.176
Autocorrelation times are between 7.65 and 11.38.
[ll/MAF:-18.338(9e+00)/ 0% | -9e-02]: 38%|███▊ | 152/400 [09:50<17:35, 4.26s/sample(s)]/home/gboehl/github/emcwrap/emcwrap/moves.py:91: RuntimeWarning: divide by zero encountered in log
np.log(sum(self.accepted)) - np.log(nchain)
[ll/MAF: 6.086(9e+00)/ 5% | -9e-02]: 50%|█████ | 201/400 [13:21<14:23, 4.34s/sample(s)]
(mcmc:) Summary from last 50 of 200 iterations (RANK, estimation tutorial):
distribution pst_mean sd/df mean sd mode hpd_5 hpd_95 \
theta beta 0.500 0.100 0.804 0.134 0.918 0.613 0.984
sigma normal 1.500 0.375 2.131 0.602 0.988 1.256 3.154
phi_pi normal 1.500 0.250 2.450 0.649 2.311 1.345 3.366
phi_y normal 0.125 0.050 0.153 0.102 0.139 -0.013 0.275
rho_u beta 0.500 0.200 0.935 0.020 0.938 0.904 0.962
rho_r beta 0.500 0.200 0.592 0.176 0.451 0.353 0.857
rho_z beta 0.500 0.200 0.788 0.345 0.458 0.052 0.999
rho beta 0.750 0.100 0.855 0.049 0.896 0.786 0.927
sig_u inv_gamma_dynare 0.100 2.000 0.238 0.094 0.201 0.108 0.385
sig_r inv_gamma_dynare 0.100 2.000 0.119 0.041 0.102 0.069 0.176
sig_z inv_gamma_dynare 0.100 2.000 0.190 0.163 0.151 0.021 0.370
error
theta 0.017
sigma 0.067
phi_pi 0.082
phi_y 0.013
rho_u 0.002
rho_r 0.022
rho_z 0.045
rho 0.006
sig_u 0.011
sig_r 0.005
sig_z 0.016
Mean acceptance fraction: 0.15
Autocorrelation times are between 9.8 and 15.2.
[ll/MAF: 8.158(9e+00)/ 0% | -3e-01]: 57%|█████▋ | 227/400 [15:11<12:18, 4.27s/sample(s)]/home/gboehl/github/emcwrap/emcwrap/moves.py:91: RuntimeWarning: divide by zero encountered in log
np.log(sum(self.accepted)) - np.log(nchain)
[ll/MAF: 10.731(9e+00)/ 5% | -1e-01]: 63%|██████▎ | 251/400 [16:53<10:25, 4.20s/sample(s)]
(mcmc:) Summary from last 50 of 250 iterations (RANK, estimation tutorial):
distribution pst_mean sd/df mean sd mode hpd_5 hpd_95 \
theta beta 0.500 0.100 0.792 0.115 0.822 0.648 0.981
sigma normal 1.500 0.375 2.286 0.540 2.554 1.469 3.195
phi_pi normal 1.500 0.250 2.304 0.483 2.579 1.270 2.825
phi_y normal 0.125 0.050 0.144 0.087 0.134 0.020 0.270
rho_u beta 0.500 0.200 0.944 0.020 0.959 0.909 0.969
rho_r beta 0.500 0.200 0.557 0.158 0.644 0.287 0.796
rho_z beta 0.500 0.200 0.842 0.318 0.981 0.169 1.000
rho beta 0.750 0.100 0.839 0.050 0.810 0.783 0.926
sig_u inv_gamma_dynare 0.100 2.000 0.212 0.108 0.168 0.098 0.431
sig_r inv_gamma_dynare 0.100 2.000 0.106 0.025 0.117 0.067 0.138
sig_z inv_gamma_dynare 0.100 2.000 0.174 0.098 0.112 0.035 0.291
error
theta 0.015
sigma 0.062
phi_pi 0.063
phi_y 0.012
rho_u 0.003
rho_r 0.019
rho_z 0.043
rho 0.006
sig_u 0.014
sig_r 0.003
sig_z 0.012
Mean acceptance fraction: 0.133
Autocorrelation times are between 11.27 and 18.85.
[ll/MAF: 10.731(9e+00)/ 0% | -5e-02]: 63%|██████▎ | 252/400 [16:57<10:19, 4.19s/sample(s)]/home/gboehl/github/emcwrap/emcwrap/moves.py:91: RuntimeWarning: divide by zero encountered in log
np.log(sum(self.accepted)) - np.log(nchain)
[ll/MAF: 15.635(5e+00)/10% | -7e-02]: 75%|███████▌ | 301/400 [20:21<06:52, 4.17s/sample(s)]
(mcmc:) Summary from last 50 of 300 iterations (RANK, estimation tutorial):
distribution pst_mean sd/df mean sd mode hpd_5 hpd_95 \
theta beta 0.500 0.100 0.787 0.081 0.973 0.668 0.954
sigma normal 1.500 0.375 2.319 0.358 2.110 1.831 2.957
phi_pi normal 1.500 0.250 2.291 0.343 1.196 1.944 2.802
phi_y normal 0.125 0.050 0.130 0.060 0.219 0.044 0.221
rho_u beta 0.500 0.200 0.951 0.015 0.957 0.928 0.971
rho_r beta 0.500 0.200 0.532 0.117 0.709 0.361 0.727
rho_z beta 0.500 0.200 0.920 0.228 0.394 0.658 0.999
rho beta 0.750 0.100 0.829 0.040 0.916 0.775 0.896
sig_u inv_gamma_dynare 0.100 2.000 0.176 0.075 0.200 0.101 0.260
sig_r inv_gamma_dynare 0.100 2.000 0.101 0.017 0.068 0.079 0.133
sig_z inv_gamma_dynare 0.100 2.000 0.152 0.075 0.131 0.044 0.254
error
theta 0.010
sigma 0.041
phi_pi 0.042
phi_y 0.007
rho_u 0.002
rho_r 0.014
rho_z 0.030
rho 0.004
sig_u 0.010
sig_r 0.002
sig_z 0.008
Mean acceptance fraction: 0.129
Autocorrelation times are between 13.0 and 22.48.
[ll/MAF: 15.212(4e+00)/ 0% | -4e-02]: 84%|████████▍ | 337/400 [22:50<04:18, 4.10s/sample(s)]/home/gboehl/github/emcwrap/emcwrap/moves.py:91: RuntimeWarning: divide by zero encountered in log
np.log(sum(self.accepted)) - np.log(nchain)
[ll/MAF: 15.067(3e+00)/15% | -2e-02]: 88%|████████▊ | 351/400 [23:47<03:22, 4.13s/sample(s)]
(mcmc:) Summary from last 50 of 350 iterations (RANK, estimation tutorial):
distribution pst_mean sd/df mean sd mode hpd_5 hpd_95 \
theta beta 0.500 0.100 0.788 0.051 0.773 0.724 0.863
sigma normal 1.500 0.375 2.233 0.287 1.784 1.772 2.633
phi_pi normal 1.500 0.250 2.301 0.277 2.489 1.982 2.716
phi_y normal 0.125 0.050 0.120 0.040 0.080 0.041 0.162
rho_u beta 0.500 0.200 0.956 0.009 0.942 0.943 0.970
rho_r beta 0.500 0.200 0.540 0.089 0.549 0.421 0.725
rho_z beta 0.500 0.200 0.977 0.113 0.998 0.990 0.999
rho beta 0.750 0.100 0.820 0.037 0.824 0.764 0.863
sig_u inv_gamma_dynare 0.100 2.000 0.154 0.035 0.187 0.097 0.196
sig_r inv_gamma_dynare 0.100 2.000 0.100 0.013 0.092 0.077 0.119
sig_z inv_gamma_dynare 0.100 2.000 0.123 0.049 0.121 0.048 0.203
error
theta 0.006
sigma 0.034
phi_pi 0.035
phi_y 0.005
rho_u 0.001
rho_r 0.011
rho_z 0.015
rho 0.005
sig_u 0.004
sig_r 0.002
sig_z 0.006
Mean acceptance fraction: 0.127
Autocorrelation times are between 14.8 and 26.49.
[ll/MAF: 15.209(2e+00)/10% | -1e-02]: 100%|██████████| 400/400 [27:09<00:00, 4.07s/sample(s)]
[save_meta:] Metadata saved as '/tmp/output/rank_tutorial_meta'
Great. On my 8 core machine this took 27 Minutes.
So where are our estimates? Our (hopefully) converged MCMC samples are currently stored in the file named rank_tutorial_sampler.h5 created by mcmc which is stored in the output_path.
You can load and use this data using the methods introduced in the processing estimation results tutorial.
[12]:
# Just for the tutorial: Cleaning the temporary directory
shutil.rmtree(output_path)