Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • giovanma/dynare
  • giorgiomas/dynare
  • Vermandel/dynare
  • Dynare/dynare
  • normann/dynare
  • MichelJuillard/dynare
  • wmutschl/dynare
  • FerhatMihoubi/dynare
  • sebastien/dynare
  • lnsongxf/dynare
  • rattoma/dynare
  • CIMERS/dynare
  • FredericKarame/dynare
  • SumuduK/dynare
  • MinjeJeon/dynare
  • camilomrch/dynare
  • DoraK/dynare
  • avtishin/dynare
  • selma/dynare
  • claudio_olguin/dynare
  • jeffjiang07/dynare
  • EthanSystem/dynare
  • stepan-a/dynare
  • wjgatt/dynare
  • JohannesPfeifer/dynare
  • gboehl/dynare
26 results
Show changes
Hide whitespace changes
Inline Side-by-side
Showing
with 745 additions and 341 deletions
examples/rbc_irf_matching.mod 0 → 100644
View file @ 9832f012
/*
* This file presents a baseline RBC model with government spending shocks where the persistence of the
* government spending shock is estimated via impulse response function (IRF) matching.
*
* Notes:
* - The empirical IRFs were estimated using the Blanchard/Perotti (2002) approach, see
* https://github.com/JohannesPfeifer/DSGE_mod/blob/master/RBC_IRF_matching/get_empirical_IRFs.m
* - They are given in the csv file rbc_irf_matching_data.csv, the first two columns contain
* the empirical IRFs of G and Y, while the third and fourth column contain the corresponding
* variances of the IRFs from a bootstrap approach.
* Importantly: this mod file does not show how to get the empirical IRFs from a SVAR model,
* but takes these as "data".
* - Of course the RBC model is not capable of generating the consumption increase
* after a government spending shock. For that reason, this mod-file only targets the IRFs for G and Y.
* - The weighting matrix uses a diagonal matrix with the inverse of the pointwise IRF variances on the main diagonal.
* - The empirical IRFs and model IRFs use an impulse size of 1 percent. Thus, there is no uncertainty about the
* initial impact. The IRF matching therefore only targets the G-response starting in the second period.
* - Note that for the current model, the number of IRFs exceeds the number of VAR parameters. Therefore,
* the distribution of the estimator will be non-standard, see Guerron-Quintana/Inoue/Kilian (2016),
* http://dx.doi.org/10.1016/j.jeconom.2016.09.009
* - The mod-file also shows how to estimate an AR(2)-process by working with the roots of the autoregressive
* process instead of the coefficients. This allows for easily restricting the process to the stability region and
* would allow specifying e.g. a beta prior for both roots as was done in Born/Peter/Pfeifer (2013), Fiscal news
* and macroeconomic volatility, https://doi.org/10.1016/j.jedc.2013.06.011
*
* Please note that the following copyright notice only applies to this Dynare
* implementation of the model.
*/
/*
* Copyright (C) 2016-17 Johannes Pfeifer,
* Copyright (C) 2024 Dynare Team
*
* This file is part of Dynare.
*
* Dynare is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Dynare is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Dynare. If not, see <https://www.gnu.org/licenses/>.
*/
%----------------------------------------------------------------
% define variables
%----------------------------------------------------------------
@#define IRF_periods=80
var y ${y}$ (long_name='output')
c ${c}$ (long_name='consumption')
k ${k}$ (long_name='capital')
l ${l}$ (long_name='hours')
z ${z}$ (long_name='TFP')
ghat ${\hat g}$ (long_name='government spending')
r ${r}$ (long_name='annualized interest rate')
w ${w}$ (long_name='real wage')
invest ${i}$ (long_name='investment')
log_y ${\log(y)}$ (long_name='log output')
log_k ${\log(k)}$ (long_name='log capital stock')
log_c ${\log(c)}$ (long_name='log consumption')
log_l ${\log(l)}$ (long_name='log labor')
log_w ${\log(w)}$ (long_name='log real wage')
log_invest ${\log(i)}$ (long_name='log investment')
;
varexo eps_z ${\varepsilon_z}$ (long_name='TFP shock')
eps_g ${\varepsilon_g}$ (long_name='government spending shock')
;
%----------------------------------------------------------------
% define parameters
%----------------------------------------------------------------
parameters
beta ${\beta}$ (long_name='discount factor')
psi ${\psi}$ (long_name='labor disutility parameter')
sigma ${\sigma}$ (long_name='risk aversion')
delta ${\delta}$ (long_name='depreciation rate')
alpha ${\alpha}$ (long_name='capital share')
rhoz ${\rho_z}$ (long_name='persistence TFP shock')
root_g_1 ${\rho_g}$ (long_name='first root of AR(2) G process')
root_g_2 ${\rho_g}$ (long_name='second root of AR(2) G process')
gammax ${\gamma_x}$ (long_name='composite growth rate')
gshare ${\frac{G}{Y}}$ (long_name='government spending share')
n ${n}$ (long_name='population growth')
x ${x}$ (long_name='technology growth (per capita output growth)')
i_y ${\frac{I}{Y}}$ (long_name='investment-output ratio')
k_y ${\frac{K}{Y}}$ (long_name='capital-output ratio')
g_ss ${\bar G}$ (long_name='government spending in steady state')
l_ss ${\bar L}$ (long_name='labor in steady state')
;
%----------------------------------------------------------------
% model equations
%----------------------------------------------------------------
model;
# rho_g_1= (root_g_1+root_g_2);
# rho_g_2= - root_g_1*root_g_2;
[name='Euler equation']
c^(-sigma)=beta/gammax*c(+1)^(-sigma)*
(alpha*exp(z(+1))*(k/l(+1))^(alpha-1)+(1-delta));
[name='Labor FOC']
psi*c^sigma*1/(1-l)=w;
[name='Law of motion capital']
gammax*k=(1-delta)*k(-1)+invest;
[name='resource constraint']
y=invest+c+g_ss*exp(ghat);
[name='production function']
y=exp(z)*k(-1)^alpha*l^(1-alpha);
[name='real wage/firm FOC labor']
w=(1-alpha)*y/l;
[name='annualized real interst rate/firm FOC capital']
r=4*alpha*y/k(-1);
[name='exogenous TFP process']
z=rhoz*z(-1)+eps_z;
[name='government spending process']
ghat=rho_g_1*ghat(-1)+rho_g_2*ghat(-2)+eps_g;
[name='Definition log output']
log_y = log(y);
[name='Definition log capital']
log_k = log(k);
[name='Definition log consumption']
log_c = log(c);
[name='Definition log hours']
log_l = log(l);
[name='Definition log wage']
log_w = log(w);
[name='Definition log investment']
log_invest = log(invest);
end;
%----------------------------------------------------------------
% set steady state values
%---------------------------------------------------------------
steady_state_model;
gammax = (1+n)*(1+x);
delta = i_y/k_y-x-n-n*x;
beta = (1+x)*(1+n)/(alpha/k_y+(1-delta));
l = l_ss;
k = ((1/beta*(1+n)*(1+x)-(1-delta))/alpha)^(1/(alpha-1))*l;
invest = (x+n+delta+n*x)*k;
y = k^alpha*l^(1-alpha);
g = gshare*y;
g_ss = g;
c = (1-gshare)*k^(alpha)*l^(1-alpha)-invest;
psi = (1-alpha)*(k/l)^alpha*(1-l)/c^sigma;
w = (1-alpha)*y/l;
r = 4*alpha*y/k;
log_y = log(y);
log_k = log(k);
log_c = log(c);
log_l = log(l);
log_w = log(w);
log_invest = log(invest);
z = 0;
ghat =0;
end;
%----------------------------------------------------------------
% calibration
%----------------------------------------------------------------
sigma = 1;
alpha = 0.33;
i_y = 0.25;
k_y = 10.4;
x = 0.0055;
n = 0.0027;
rhoz = 0.97;
root_g_1 = 0.9602;
root_g_2 = 0;
gshare = 0.2038;
l_ss = 1/3;
shocks;
var eps_g = 1;
end;
steady;
check;
varobs ghat log_y y; // you need to specify observables
%----------------------------------------------------------------
% IRF matching example 1:
% - different ways to MANUALLY enter values and weights
% - Maximum likelihood estimation
%----------------------------------------------------------------
estimated_params;
root_g_1 , 0.90 , 0, 1;
root_g_2 , 0.10 , 0, 1;
end;
xx = [1.007,1.117,1.092];
ww = [51,52];
matched_irfs;
var log_y ; varexo eps_g ; periods 1, 2 ; values 0.20, 0.17 ; weights 360, 140 ;
var ghat ; varexo eps_g ; periods 2 3:5 ; values 1.01, (xx) ; weights 50, 20 ;
var y ; varexo eps_g ; periods 10:11 ; values (log(1.05)) ; weights (ww) ;
end;
method_of_moments(mom_method = irf_matching, mode_compute = 5, additional_optimizer_steps=[4]);
%----------------------------------------------------------------
% IRF matching example 2
% - use all IRFs given in MATLAB objects
% - use Bayesian Slice sampler
%----------------------------------------------------------------
estimated_params(overwrite);
root_g_1 , 0.50 , 0, 1, beta_pdf , 0.50 , 0.20;
root_g_2 , 0.10 , 0, 1, beta_pdf , 0.50 , 0.20;
end;
% get data
irfs_data = importdata('rbc_irf_matching_data.csv');
irfs_ghat_eps_g = irfs_data.data(2:80,1); % start in t=2 due to identification restrictions in SVAR
irfs_log_y_eps_g = irfs_data.data(1:80,2);
weights_ghat_eps_g = 1./irfs_data.data(2:80,3);
weights_log_y_eps_g = 1./irfs_data.data(1:80,4);
matched_irfs(overwrite);
var ghat ; varexo eps_g; periods 2:80; values (irfs_ghat_eps_g); weights (weights_ghat_eps_g);
var log_y; varexo eps_g; periods 1:80; values (irfs_log_y_eps_g); weights (weights_log_y_eps_g);
end;
method_of_moments(mom_method = irf_matching
,order = 1
,mh_nblocks = 2, mh_replic = 50
,posterior_sampling_method = 'slice'
,plot_priors = 1
);
%----------------------------------------------------------------
% IRF matching example 3:
% - use anonymous function to access IRFs more flexibly
% - showcase how to use irf_matching_file
% - find posterior mode
%----------------------------------------------------------------
% get data
irfs_data = importdata('rbc_irf_matching_data.csv');
% use anonymous function (or MATLAB function) to have more flexibility, but inputs can only be numerical
% we also take 100 just for illustration that you can do any required transformation in an irf_matching_file
irfs_vals = @(j) 100.*(irfs_data.data(2:80,j));
irfs_weights = @(j) 1./(irfs_data.data(2:80,j));
matched_irfs(overwrite);
var ghat ; varexo eps_g; periods 2:80; values (irfs_vals(1)); weights (irfs_weights(3));
var log_y; varexo eps_g; periods 2:80; values (irfs_vals(2)); weights (irfs_weights(4));
end;
% we use the irf_matching_file to transform variable y to log(y) so the model
% variable aligns with the variable from the given empirical SVAR
method_of_moments(mom_method = irf_matching
,irf_matching_file = rbc_irf_matching_transformations
,mh_replic = 0,plot_priors = 0
);
\ No newline at end of file
examples/rbc_irf_matching_data.csv 0 → 100644
View file @ 9832f012
G,Y,var_G,var_Y
1,0.208070478344545,1.12054105855257e-33,0.0027602236508546
1.01463693816715,0.166450236530301,0.00583152879723266,0.00679133276154495
1.00764238152471,0.230433660500946,0.0111184998994349,0.00952654639548326
1.11696510115769,0.215252336067127,0.0164663222282161,0.0113464752654877
1.09223998892147,0.222958847463141,0.0216966728235886,0.013533420700953
1.03736643322521,0.243972174559692,0.0257004687001329,0.0151053406229588
1.01711112212948,0.261016995864984,0.0303544421773506,0.0162405771811019
0.967577038730636,0.269060130173636,0.0339076793971577,0.0171352445593892
0.912649802559854,0.278584624341349,0.03620212016322,0.0179592573395723
0.867058516612021,0.283001682449059,0.038151205713992,0.0187783728830758
0.820597806619289,0.283598412003388,0.0396511019185971,0.0194594651103199
0.77563492323912,0.281760898055785,0.0403347277351183,0.0198713886600475
0.735846344956036,0.277110769410598,0.0407036832985301,0.0198800923047385
0.698850036380633,0.269931336693529,0.0407917065478011,0.0194569744071588
0.664596916513962,0.261143278620547,0.0405012612908696,0.0186271529903141
0.633528960929983,0.251064751428521,0.0400022455114944,0.0174559567202935
0.604912703413876,0.240132445806634,0.0393594337883548,0.0160603819816583
0.578360605178466,0.228823340604282,0.0385445703756185,0.0145721875285327
0.55371549942988,0.217458827167635,0.0376075435183867,0.0130962012473614
0.530621264922238,0.206298301306342,0.0365729915386994,0.0117094186557471
0.508813730835084,0.195562831928213,0.0354403165776896,0.010461528435305
0.488125279535286,0.185390333897221,0.0342268285355491,0.00937102041748079
0.468386758403429,0.175863201796663,0.0329491563166573,0.00843340667968603
0.449471041462416,0.167022878990409,0.0316223984500399,0.00763260719022832
0.431296406103216,0.158870322365946,0.0302644814355948,0.00694575471751527
0.413799596231625,0.151377829877172,0.0288934660918837,0.00634988691182681
0.396937751823371,0.144499570770831,0.0275270293690654,0.00582564616104379
0.380685985677042,0.13817717429653,0.0261806418593218,0.00535836547483349
0.365028340579903,0.132346351995696,0.0248665732796567,0.00493799567463183
0.349954881191496,0.126942318129237,0.0235951459469029,0.0045584177486503
0.335459383987806,0.121903145354258,0.022373467574597,0.00421578980579863
0.321536118675061,0.117172332783573,0.021205727363848,0.00390754150917342
0.308178316488601,0.112700463782387,0.0200943664605458,0.00363147379277453
0.295377351869768,0.108445858384018,0.019040123479836,0.00338513596160304
0.283122066100135,0.104374632848495,0.0180423053448312,0.00316565728079348
0.271398650164465,0.100460305211796,0.0170995738403303,0.0029698495549103
0.260190863606074,0.0966830386847356,0.0162100638803249,0.00279435029276752
0.249480370591311,0.0930286987898092,0.0153715252985084,0.00263589300567087
0.239247213436084,0.0894878338463341,0.0145815957210218,0.002491545439791
0.229470346306504,0.0860546497623421,0.0138378615720872,0.00235884935269597
0.220128148976686,0.0827260491927552,0.0131378639605267,0.00223588384965084
0.211198902955524,0.0795007779016462,0.0124791808709379,0.00212125323555392
0.20266120682353,0.0763786992809964,0.0118594333328108,0.0020139945677632
0.194494311443246,0.0733602078340912,0.0112762767066973,0.00191345670528973
0.186678374172779,0.0704457804639186,0.0107274179119509,0.0018191767154505
0.179194635081175,0.0676356552603192,0.0102106211287206,0.00173077210770515
0.172025520775263,0.064929622999291,0.00972370794405955,0.0016478667536602
0.165154686254735,0.0623269135604257,0.00926457432688479,0.00157005669764952
0.158567006362608,0.0598261581740075,0.00883120737460965,0.00149690646906453
0.152248528236286,0.0574254089705774,0.00842170119327602,0.00142796690804853
0.146186396016794,0.0551221987689224,0.00803427530290723,0.0013628024472437
0.140368757967879,0.0529136261066764,0.00766728912012237,0.00130101642649392
0.134784664597338,0.0507964530245216,0.00731924890008266,0.00124226770384004
0.129423964825334,0.0487672057100314,0.00698880917814805,0.00118627695274624
0.124277205618631,0.046822270602385,0.00667476765058639,0.00113282353216918
0.119335538941298,0.0449579808540134,0.0063760543634253,0.0010817364372413
0.11459063850232,0.0431706900378893,0.00609171803409054,0.00103288330782315
0.110034627608638,0.0414568316344016,0.00582091131964085,0.000986160390158817
0.105660018476945,0.0398129641266707,0.00556287623442419,0.000941484872060258
0.101459662633845,0.0382358024894679,0.00531693108962639,0.000898789756847698
0.0974267115198392,0.0367222375053657,0.00508245926485931,0.000858020371449084
0.0935545860830807,0.0352693447261226,0.00485889955647229,0.000819131356725285
0.0898369539782651,0.0338743850644921,0.00464573784815541,0.00078208329876107
0.0862677129419094,0.0325347989994978,0.00444249975285543,0.000746838714885712
0.0828409789666302,0.0312481962526895,0.00424874389794858,0.000713357679487274
0.0795510780159448,0.0300123425850257,0.00406405581698042,0.000681593780987298
0.076392540182423,0.0288251451088667,0.00388804253930626,0.000651491166504331
0.0733600953745436,0.0276846372357925,0.00372032800246882,0.000622983230006533
0.070448669804769,0.0265889641106603,0.00356054943074244,0.000595993138422716
0.0676533827305797,0.0255363691313685,0.0034083547518209,0.00057043599416774
0.064969543062731,0.024525181932661,0.00326340099829315,0.00054622212266834
0.0623926455954874,0.0235538080267355,0.00312535355830137,0.000523260827093126
0.0599183667294525,0.0226207201454146,0.00299388608378455,0.000501463968548093
0.0575425596484649,0.0217244512173301,0.00286868083997892,0.000480748874949654
0.055261248979131,0.0208635888360229,0.00274942930127684,0.000461040295246497
0.0530706250072852,0.0200367710268393,0.00263583284298659,0.000442271333040507
0.0509670375530574,0.0192426830970757,0.00252760342283966,0.000424383467067322
0.0489469896186406,0.0184800553498327,0.00242446418777184,0.000407325873359123
0.0470071309236527,0.0177476614525368,0.0023261499716701,0.000391054302991331
0.0451442514353558,0.0170443172715253,0.00223240766511842,0.000375529754828271
examples/rbc_irf_matching_transformations.m 0 → 100644
View file @ 9832f012
function [modelIrf, error_indicator] = rbc_irf_matching_transformations(modelIrf, M_, options_mom_, ys_)
% -------------------------------------------------------------------------
% This file manipulates model IRFs to be consistent with empirical IRFS
% -------------------------------------------------------------------------
% INPUTS
% - modelIrf: [options_mom_.irf by M_.endo_nbr by M_.exo_nbr]
% array of IRFs for all model variables and all shocks
% - M_: [structure] Dynare model structure
% - options_mom_: [structure] Dynare options structure
% - ys_: [double] steady state values of all endogenous variables
% -------------------------------------------------------------------------
% OUTPUTS
% - modelIrf: [options_mom_.irf by M_.endo_nbr by M_.exo_nbr]
% modified array of IRFs for all model variables and all shocks
% - error_indicator: [boolean] indicator of success (0) or failure (1)
% -------------------------------------------------------------------------
% This function is called by
% - mom.run
% -------------------------------------------------------------------------
% Copyright © 2024 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
% initialize error indicator
error_indicator = 0;
% get indices of variables
idx_ghat = find(ismember(M_.endo_names,'ghat'));
idx_log_y = find(ismember(M_.endo_names,'log_y'));
idx_eps_g = find(ismember(M_.exo_names,'eps_g'));
% manipulate the model IRFs to match the empirical IRFs (e.g. cumsum, common scaling, trends, ratios, etc.)
modelIrf(:,idx_ghat,idx_eps_g) = 100.*modelIrf(:,idx_ghat,idx_eps_g);
modelIrf(:,idx_log_y,idx_eps_g) = 100.*modelIrf(:,idx_log_y,idx_eps_g);
end
\ No newline at end of file
license.txt
View file @ 9832f012
Format: https://www.debian.org/doc/packaging-manuals/copyright-format/1.0/
Upstream-Name: Dynare
Upstream-Contact: Dynare Team, whose members in 2023 are:
Upstream-Contact: Dynare Team, whose members in 2025 are:
- Stéphane Adjemian <stephane.adjemian@univ-lemans.fr>
- Michel Juillard <michel.juillard@mjui.fr>
- Frédéric Karamé <frederic.karame@univ-lemans.fr>
......@@ -23,7 +23,7 @@ Upstream-Contact: Dynare Team, whose members in 2023 are:
Source: https://www.dynare.org
Files: *
Copyright: 1996-2023 Dynare Team
Copyright: 1996-2025 Dynare Team
License: GPL-3+
Files: matlab/+occbin/IVF_core.m
......@@ -86,13 +86,13 @@ License: public-domain-aim
Journal of Economic Dynamics and Control, 2010, vol. 34, issue 3,
pages 472-489
Files: matlab/optimization/bfgsi1.m matlab/csolve.m matlab/optimization/csminit1.m matlab/optimization/numgrad2.m
Files: matlab/optimization/bfgsi1.m matlab/optimization/csolve.m matlab/optimization/csminit1.m matlab/optimization/numgrad2.m
matlab/optimization/numgrad3.m matlab/optimization/numgrad3_.m matlab/optimization/numgrad5.m
matlab/optimization/numgrad5_.m matlab/optimization/csminwel1.m matlab/+bvar/density.m
matlab/+bvar/toolbox.m matlab/partial_information/PI_gensys.m matlab/partial_information/qzswitch.m
matlab/partial_information/qzdiv.m
matlab/+bvar/toolbox.m matlab/partial_information/PI_gensys.m matlab/partial_information/PI_qzswitch.m
matlab/partial_information/PI_qzdiv.m
Copyright: 1993-2009 Christopher Sims
2006-2023 Dynare Team
2006-2024 Dynare Team
License: GPL-3+
Files: matlab/optimization/cmaes.m
......@@ -123,6 +123,11 @@ Copyright: 2010-2015 Alexander Meyer-Gohde
2015-2017 Dynare Team
License: GPL-3+
Files: matlab/collapse_figures_in_tabgroup.m
Copyright: 2023 Eduard Benet Cerda
2024 Dynare Team
License: GPL-3+
Files: matlab/convergence_diagnostics/raftery_lewis.m
Copyright: 2016 Benjamin Born and Johannes Pfeifer
2016-2017 Dynare Team
......@@ -172,7 +177,7 @@ Comment: Written by Jessica Cariboni and Francesca Campolongo
Files: matlab/+gsa/cumplot.m
matlab/+gsa/monte_carlo_filtering.m
matlab/+gsa/skewness.m
matlab/+gsa/log_trans_.m
matlab/+gsa/log_transform.m
matlab/+gsa/map_calibration.m
matlab/+gsa/map_identification.m
matlab/+gsa/monte_carlo_filtering_analysis.m
......@@ -247,11 +252,11 @@ License: BSD-2-clause
Files: examples/fs2000_data.m
Copyright: 2000-2022 Frank Schorfheide
Copyright: 2023 Dynare Team
2023 Dynare Team
License: CC-BY-SA-4.0
Files: doc/*.rst doc/*.tex doc/*.svg doc/*.pdf doc/*.bib
Copyright: 1996-2022 Dynare Team
Copyright: 1996-2025 Dynare Team
License: GFDL-NIV-1.3+
Files: doc/dr.tex doc/dr.bib
......@@ -292,28 +297,6 @@ Files: preprocessor/doc/preprocessor/*
Copyright: 2007-2019 Dynare Team
License: CC-BY-SA-4.0
Files: contrib/jsonlab/*
Copyright: 2011-2020 Qianqian Fang <q.fang at neu.edu>
2016 Bastian Bechtold
License: GPL-3+ or BSD-3-clause
Files: contrib/jsonlab/base64decode.m
contrib/jsonlab/base64encode.m
contrib/jsonlab/gzipdecode.m
contrib/jsonlab/gzipencode.m
contrib/jsonlab/zlibdecode.m
contrib/jsonlab/zlibencode.m
Copyright: 2012 Kota Yamaguchi
2011-2020 Qianqian Fang <q.fang at neu.edu>
License: GPL-3+ or BSD-2-clause
Files: contrib/jsonlab/loadjson.m
Copyright: 2011-2020 Qianqian Fang
2009 Nedialko Krouchev
2009 François Glineur
2008 Joel Feenstra
License: GPL-3+ or BSD-2-clause or BSD-3-clause
Files: contrib/ms-sbvar/utilities_dw/*
Copyright: 1996-2011 Daniel Waggoner
License: GPL-3+
......@@ -420,32 +403,6 @@ License: BSD-2-clause
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
POSSIBILITY OF SUCH DAMAGE.
License: BSD-3-clause
Redistribution and use in source and binary forms, with or without modification,
are permitted provided that the following conditions are met:
.
* Redistributions of source code must retain the above copyright notice, this
list of conditions and the following disclaimer.
.
* Redistributions in binary form must reproduce the above copyright notice,
this list of conditions and the following disclaimer in the documentation
and/or other materials provided with the distribution.
.
* Neither the name of the copyright holder nor the names of its contributors
may be used to endorse or promote products derived from this software without
specific prior written permission.
.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
License: GFDL-NIV-1.3+
Permission is granted to copy, distribute and/or modify this document
under the terms of the GNU Free Documentation License, Version 1.3 or
......
macOS/build.sh
View file @ 9832f012
#!/usr/bin/env bash
# Copyright © 2019-2023 Dynare Team
# Copyright © 2019-2024 Dynare Team
#
# This file is part of Dynare.
#
......@@ -34,11 +34,16 @@ if [[ "$PKG_ARCH" == arm64 ]]; then
MATLAB_ARCH=maca64
else
BREWDIR=/usr/local
# Remove /opt/homebrew/bin from PATH, so it does not intervene with the x86_64 compilations
# Remove /opt/homebrew/bin from PATH, so it does not interfere with the x86_64 compilations
path_remove PATH /opt/homebrew/bin
MATLAB_ARCH=maci64
fi
MATLAB_PATH=/Applications/"$PKG_ARCH"/MATLAB_R2023b.app
MATLAB_PATH=/Applications/"$PKG_ARCH"/MATLAB_R2024b.app
# Workaround for bug in Xcode 15.3 which does not include m4
# See https://github.com/Homebrew/homebrew-core/issues/165388
# and https://trac.macports.org/ticket/69639
path_prepend PATH "$BREWDIR"/opt/m4/bin
# Append texbin to PATH to access latexmk and friends
path_prepend PATH /Library/TeX/texbin
......@@ -65,18 +70,18 @@ ln -s "$BREWDIR"/opt/gcc/lib/gcc/"$GCC_VERSION"/libquadmath.a "$QUADMATH_DIR"
cd "$ROOTDIR"
# NB: the addition of -Wl,-ld_classic is a workaround for https://github.com/mesonbuild/meson/issues/12282 (see also the native file)
common_meson_opts=(-Dbuild_for=matlab -Dbuildtype=release -Dprefer_static=true -Dfortran_args="[ '-B', '$LIB64/Slicot/' ]" \
common_meson_opts=(-Dbuild_for=matlab --buildtype=release --prefer-static -Dfortran_args="[ '-B', '$LIB64/Slicot/' ]" \
-Dc_link_args="[ '-Wl,-ld_classic', '-L$QUADMATH_DIR' ]" -Dcpp_link_args="[ '-Wl,-ld_classic', '-L$QUADMATH_DIR' ]" -Dfortran_link_args="[ '-Wl,-ld_classic', '-L$QUADMATH_DIR' ]" \
--native-file macOS/homebrew-native-$PKG_ARCH.ini)
# Build for MATLAB ⩾ R2018b (x86_64) and MATLAB ⩾ R2023b (arm64)
arch -"$PKG_ARCH" meson setup "${common_meson_opts[@]}" -Dmatlab_path="$MATLAB_PATH" build-matlab --wipe
arch -"$PKG_ARCH" meson compile -v -C build-matlab
arch -"$PKG_ARCH" meson setup "${common_meson_opts[@]}" -Dmatlab_path="$MATLAB_PATH" build-macOS-matlab --wipe
arch -"$PKG_ARCH" meson compile -v -C build-macOS-matlab
# If not in CI, build the docs
if [[ -z $CI ]]; then
arch -"$PKG_ARCH" meson compile -v -C build-matlab doc
ln -s build-matlab build-doc
arch -"$PKG_ARCH" meson compile -v -C build-macOS-matlab doc
ln -s build-macOS-matlab build-doc
fi
##
......@@ -85,7 +90,7 @@ fi
# Determine Dynare version if not passed by an environment variable as in the CI
if [[ -z $VERSION ]]; then
cd build-matlab
cd build-macOS-matlab
VERSION=$(meson introspect --projectinfo | sed -En 's/^.*"version": "([^"]*)".*$/\1/p')
cd ..
fi
......@@ -117,9 +122,9 @@ mkdir -p \
"$PKGFILES"/scripts \
"$PKGFILES"/contrib/ms-sbvar/TZcode
if [[ "$PKG_ARCH" == x86_64 ]]; then
mkdir -p "$PKGFILES"/mex/matlab/"$MATLAB_ARCH"-9.5-23.2
mkdir -p "$PKGFILES"/mex/matlab/"$MATLAB_ARCH"-9.5-24.2
else
mkdir -p "$PKGFILES"/mex/matlab/"$MATLAB_ARCH"-23.2
mkdir -p "$PKGFILES"/mex/matlab/"$MATLAB_ARCH"-23.2-24.2
fi
cp -p "$ROOTDIR"/NEWS.md "$PKGFILES"
......@@ -127,20 +132,20 @@ cp -p "$ROOTDIR"/COPYING "$PKGFILES"
cp -p "$ROOTDIR"/license.txt "$PKGFILES"
cp -pr "$ROOTDIR"/matlab "$PKGFILES"
cp -p "$ROOTDIR"/build-matlab/dynare_version.m "$PKGFILES"/matlab
cp -p "$ROOTDIR"/build-macOS-matlab/dynare_version.m "$PKGFILES"/matlab
cp -pr "$ROOTDIR"/examples "$PKGFILES"
cp -p "$ROOTDIR"/build-matlab/preprocessor/src/dynare-preprocessor "$PKGFILES"/preprocessor
cp -p "$ROOTDIR"/build-macOS-matlab/preprocessor/src/dynare-preprocessor "$PKGFILES"/preprocessor
# Create backward-compatibility symlink
mkdir -p "$PKGFILES"/matlab/preprocessor64
ln -sf ../../preprocessor/dynare-preprocessor "$PKGFILES"/matlab/preprocessor64/dynare_m
if [[ "$PKG_ARCH" == x86_64 ]]; then
cp -L "$ROOTDIR"/build-matlab/*.mex"$MATLAB_ARCH" "$PKGFILES"/mex/matlab/"$MATLAB_ARCH"-9.5-23.2
cp -L "$ROOTDIR"/build-macOS-matlab/*.mex"$MATLAB_ARCH" "$PKGFILES"/mex/matlab/"$MATLAB_ARCH"-9.5-24.2
else
cp -L "$ROOTDIR"/build-matlab/*.mex"$MATLAB_ARCH" "$PKGFILES"/mex/matlab/"$MATLAB_ARCH"-23.2
cp -L "$ROOTDIR"/build-macOS-matlab/*.mex"$MATLAB_ARCH" "$PKGFILES"/mex/matlab/"$MATLAB_ARCH"-23.2-24.2
fi
cp -p "$ROOTDIR"/scripts/dynare.el "$PKGFILES"/scripts
......@@ -150,8 +155,8 @@ cp "$ROOTDIR"/build-doc/*.pdf "$PKGFILES"
cp "$ROOTDIR"/build-doc/preprocessor/doc/*.pdf "$PKGFILES"/doc
cp -r "$ROOTDIR"/build-doc/dynare-manual.html "$PKGFILES"/doc
mkdir -p "$PKGFILES"/matlab/modules/dseries/externals/x13/macOS/64
cp -p "$ROOTDIR"/macOS/deps/"$PKG_ARCH"/lib64/x13as/x13as "$PKGFILES"/matlab/modules/dseries/externals/x13/macOS/64
mkdir -p "$PKGFILES"/matlab/dseries/externals/x13/macOS/64
cp -p "$ROOTDIR"/macOS/deps/"$PKG_ARCH"/lib64/x13as/x13as "$PKGFILES"/matlab/dseries/externals/x13/macOS/64
cd "$ROOTDIR"/macOS/pkg
......
macOS/homebrew-native-arm64.ini
View file @ 9832f012
# Meson native file for compiling under Homebrew for arm64 architecture
[binaries]
cpp = '/opt/homebrew/bin/g++-13'
c = '/opt/homebrew/bin/gcc-13'
cpp = '/opt/homebrew/bin/g++-14'
c = '/opt/homebrew/bin/gcc-14'
flex = '/opt/homebrew/opt/flex/bin/flex'
bison = '/opt/homebrew/opt/bison/bin/bison'
......
macOS/homebrew-native-x86_64.ini
View file @ 9832f012
# Meson native file for compiling under Homebrew for x86_64 architecture
[binaries]
cpp = '/usr/local/bin/g++-13'
c = '/usr/local/bin/gcc-13'
cpp = '/usr/local/bin/g++-14'
c = '/usr/local/bin/gcc-14'
flex = '/usr/local/opt/flex/bin/flex'
bison = '/usr/local/opt/bison/bin/bison'
......
matlab/+conditional_forecasts/run.m
View file @ 9832f012
......@@ -80,7 +80,7 @@ end
if estimated_model
if options_.prefilter
error('imcforecast:: Conditional forecasting does not support the prefiltering option')
error('conditional_forecast: Conditional forecasting does not support the prefiltering option')
end
if ischar(options_cond_fcst.parameter_set)
switch options_cond_fcst.parameter_set
......@@ -103,19 +103,19 @@ if estimated_model
xparam = bayestopt_.p1;
graph_title='Prior Mean';
otherwise
disp('imcforecast:: If the input argument is a string, then it has to be equal to:')
disp(' ''calibration'', ')
disp(' ''posterior_mode'', ')
disp(' ''posterior_mean'', ')
disp(' ''posterior_median'', ')
disp(' ''prior_mode'' or')
disp(' ''prior_mean''.')
error('imcforecast:: Wrong argument type!')
disp('conditional_forecast: If the input argument is a string, then it has to be equal to:')
disp(' ''calibration'', ')
disp(' ''posterior_mode'', ')
disp(' ''posterior_mean'', ')
disp(' ''posterior_median'', ')
disp(' ''prior_mode'' or')
disp(' ''prior_mean''.')
error('conditional_forecast: Wrong argument type!')
end
else
xparam = options_cond_fcst.parameter_set;
if length(xparam)~=length(M_.params)
error('imcforecast:: The dimension of the vector of parameters doesn''t match the number of estimated parameters!')
error('conditional_forecast: The dimension of the vector of parameters doesn''t match the number of estimated parameters!')
end
end
set_parameters(xparam);
......@@ -219,7 +219,7 @@ n2 = size(options_cond_fcst.controlled_varexo,1);
constrained_vars = oo_.dr.inv_order_var(constrained_vars); % must be in decision rule order
if n1 ~= n2
error('imcforecast:: The number of constrained variables doesn''t match the number of controlled shocks')
error('conditional_forecast: The number of constrained variables doesn''t match the number of controlled shocks')
end
% Get indices of controlled varexo.
......
matlab/+gsa/run.m
View file @ 9832f012
......@@ -13,7 +13,7 @@ function x0=run(M_,oo_,options_,bayestopt_,estim_params_,options_gsa)
% M. Ratto (2008), Analysing DSGE Models with Global Sensitivity Analysis,
% Computational Economics (2008), 31, pp. 115–139
% Copyright © 2008-2023 Dynare Team
% Copyright © 2008-2024 Dynare Team
%
% This file is part of Dynare.
%
......@@ -101,6 +101,9 @@ if ~isempty(options_gsa.datafile) || isempty(bayestopt_) || options_gsa.rmse
disp('must be specified for RMSE analysis!');
error('Sensitivity anaysis error!')
end
if isfield(options_gsa,'nobs')
options_.nobs=options_gsa.nobs;
end
if ~isempty(options_.nobs) && length(options_.nobs)~=1
error('dynare_sensitivity does not support recursive estimation. Please specify nobs as a scalar, not a vector.')
end
......@@ -108,9 +111,6 @@ if ~isempty(options_gsa.datafile) || isempty(bayestopt_) || options_gsa.rmse
if isfield(options_gsa,'first_obs')
options_.first_obs=options_gsa.first_obs;
end
if isfield(options_gsa,'nobs')
options_.nobs=options_gsa.nobs;
end
if isfield(options_gsa,'presample')
options_.presample=options_gsa.presample;
end
......@@ -150,6 +150,11 @@ end
[~,~,~,~,oo_.dr,M_.params] = dynare_resolve(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
if isfield(oo_.dr,'eigval') && any(abs(oo_.dr.eigval-1)<abs(1-options_.qz_criterium)) && options_.qz_criterium<1
fprintf('\ngsa: The model features a unit root, but qz_criterium<1. Check whether that is intended.')
fprintf('\ngsa: If not, use the diffuse_filter option.\n')
end
options_gsa = set_default_option(options_gsa,'identification',0);
if options_gsa.identification
options_gsa.redform=0;
......@@ -355,7 +360,7 @@ if options_gsa.redform && ~isempty(options_gsa.namendo)
if isempty(options_gsa.threshold_redform) && ~(exist('gsa_sdp','file')==6 || exist('gsa_sdp','file')==2)
fprintf('\nThe "SS-ANOVA-R: MATLAB Toolbox for the estimation of Smoothing Spline ANOVA models with Recursive algorithms" is missing.\n')
fprintf('To obtain it, go to:\n\n')
fprintf('https://ec.europa.eu/jrc/en/macro-econometric-statistical-software/ss-anova-r-downloads \n\n')
fprintf('https://joint-research-centre.ec.europa.eu/system/files/2025-01/ss_anova_recurs.zip \n\n')
fprintf('and follow the instructions there.\n')
fprintf('After obtaining the files, you need to unpack them and set a Matlab Path to those files.\n')
error('SS-ANOVA-R Toolbox missing!')
......@@ -522,4 +527,4 @@ if options_gsa.glue
save([OutputDirectoryName,'/',fname_,'_glue_mc.mat'], 'Out', 'Sam', 'Lik', 'Obs', 'Rem','Info', 'Exo')
end
end
end
\ No newline at end of file
end
matlab/+identification/analysis.m
View file @ 9832f012
......@@ -297,10 +297,19 @@ if info(1) == 0 %no errors in solution
options_.analytic_derivation = -2; %this sets asy_Hess=1 in dsge_likelihood.m
[info, oo_, options_, M_] = stoch_simul(M_, options_, oo_, options_.varobs);
dataset_ = dseries(oo_.endo_simul(options_.varobs_id,100+1:end)',dates('1Q1'), options_.varobs); %get information on moments
% set info on missing data
if dataset_info.missing.state
[dataset_info.missing.aindex, dataset_info.missing.number_of_observations, dataset_info.missing.no_more_missing_observations, dataset_info.missing.vindex] = ...
describe_missing_data(dataset_.data);
else
dataset_info.missing.aindex = num2cell(transpose(repmat(1:dataset_.vobs,dataset_.nobs,1)),1);
dataset_info.missing.no_more_missing_observations = 1;
end
derivatives_info.no_DLIK = 1;
bounds = prior_bounds(bayestopt_, options_.prior_trunc); %reset bounds as lb and ub must only be operational during mode-finding
%note that for order>1 we do not provide any information on DT,DYss,DOM in derivatives_info, such that dsge_likelihood creates an error. Therefore the computation will be based on simulated_moment_uncertainty for order>1.
[~, info, ~, ~, AHess, ~, ~, M_, options_, ~, oo_.dr] = dsge_likelihood(params', dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, bounds, oo_.dr, oo_.steady_state,oo_.exo_steady_state, oo_.exo_det_steady_state. derivatives_info); %non-used output variables need to be set for octave for some reason
[~, info, ~, ~, AHess, ~, ~, M_, options_, ~, oo_.dr] = dsge_likelihood(params', dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, bounds, oo_.dr, oo_.steady_state,oo_.exo_steady_state, oo_.exo_det_steady_state, derivatives_info); %non-used output variables need to be set for octave for some reason
%note that for the order of parameters in AHess we have: stderr parameters come first, corr parameters second, model parameters third. the order within these blocks corresponds to the order specified in the estimated_params block
options_.analytic_derivation = analytic_derivation; %reset option
AHess = -AHess; %take negative of hessian
......
matlab/+identification/run.m
View file @ 9832f012
......@@ -74,6 +74,7 @@ fname = M_.fname; %model name
dname = M_.dname; %model name
%turn warnings off, either globally or only relevant ids
orig_warning_state = warning;
if isoctave
%warning('off'),
warning('off','Octave:singular-matrix');
......@@ -239,7 +240,7 @@ end
% check for external draws, i.e. set pdraws0 for a gsa analysis
if options_ident.gsa_sample_file
GSAFolder = checkpath('gsa',dname);
GSAFolder = CheckPath('gsa',dname);
if options_ident.gsa_sample_file==1
load([GSAFolder,filesep,fname,'_prior'],'lpmat','lpmat0','istable');
elseif options_ident.gsa_sample_file==2
......@@ -484,7 +485,7 @@ if iload <=0
kk=0;
while kk<50 && info(1)
kk=kk+1;
params = Prior.draw();
params = Prior.draw()'; %column vector is expected
options_ident.tittxt = 'Random_prior_params'; %title text for graphs and figures
% perform identification analysis
[ide_moments_point, ide_spectrum_point, ide_minimal_point, ide_hess_point, ide_reducedform_point, ide_dynamic_point, ~, info, error_indicator_point] = ...
......@@ -498,6 +499,7 @@ if iload <=0
fprintf('The model did not solve for any of 50 attempts of random samples from the prior\n');
end
fprintf('-----------\n');
warning(orig_warning_state);
return
else
% found a (random) point that solves the model
......@@ -976,6 +978,6 @@ if SampleSize > 1
end
%reset warning state
warning_config;
warning(orig_warning_state);
fprintf('\n==== Identification analysis completed ====\n\n')
matlab/+identification/simulated_moment_uncertainty.m
View file @ 9832f012
......@@ -12,7 +12,7 @@ function [cmm, mm] = simulated_moment_uncertainty(indx, periods, replic,options_
% - cmm: [n_moments by n_moments] covariance matrix of simulated moments
% - mm: [n_moments by replic] matrix of moments
% Copyright © 2009-2018 Dynare Team
% Copyright © 2009-2024 Dynare Team
%
% This file is part of Dynare.
%
......@@ -60,15 +60,10 @@ end
oo_.dr=set_state_space(oo_.dr,M_);
if options_.logged_steady_state %if steady state was previously logged, undo this
oo_.dr.ys=exp(oo_.dr.ys);
oo_.steady_state=exp(oo_.steady_state);
options_.logged_steady_state=0;
logged_steady_state_indicator=1;
evalin('base','options_.logged_steady_state=0;')
else
logged_steady_state_indicator=0;
end
[oo_.dr,info,M_.params] = compute_decision_rules(M_,options_,oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
......@@ -92,7 +87,6 @@ else
end
end
for j=1:replic
[ys, oo_.exo_simul] = simult(y0,oo_.dr,M_,options_);%do simulation
oo_=disp_moments(ys, options_.varobs,M_,options_,oo_); %get moments
......@@ -106,8 +100,5 @@ for j=1:replic
end
dyn_waitbar_close(h);
if logged_steady_state_indicator
evalin('base','options_.logged_steady_state=1;') %reset base workspace option to conform to base oo_
end
cmm = cov(mm');
disp('Simulated moment uncertainty ... done!')
matlab/+mom/Jtest.m
View file @ 9832f012
......@@ -70,4 +70,6 @@ if options_mom_.mom.mom_nbr > length(xparam)
J_test.p_val = 1-chi2cdf(J_test.j_stat, J_test.degrees_freedom);
fprintf('\nValue of J-test statistic: %f\n',J_test.j_stat);
fprintf('p-value of J-test statistic: %f\n',J_test.p_val);
else
J_test=[];
end
\ No newline at end of file
matlab/+mom/default_option_mom_values.m
View file @ 9832f012
function options_mom_ = default_option_mom_values(options_mom_, options_, dname, do_bayesian_estimation)
% options_mom_ = default_option_mom_values(options_mom_, options_, dname, do_bayesian_estimation)
function options_mom_ = default_option_mom_values(options_mom_, options_, dname, fname, do_bayesian_estimation)
% options_mom_ = default_option_mom_values(options_mom_, options_, dname, fname, do_bayesian_estimation)
% -------------------------------------------------------------------------
% Returns structure containing the options for method_of_moments command.
% Note 1: options_mom_ is local and contains default and user-specified
......@@ -16,6 +16,7 @@ function options_mom_ = default_option_mom_values(options_mom_, options_, dname,
% o options_mom_: [structure] all user-specified settings (from the method_of_moments command)
% o options_: [structure] global options
% o dname: [string] default name of directory to store results
% o fname: [string] default name of mod file
% o do_bayesian_estimation [boolean] indicator whether we do Bayesian estimation
% -------------------------------------------------------------------------
% OUTPUTS
......@@ -177,7 +178,7 @@ options_mom_ = set_default_option(options_mom_,'dr_cycle_reduction',false);
options_mom_ = set_default_option(options_mom_,'dr_cycle_reduction_tol',1e-7); % convergence criterion used in the cycle reduction algorithm
options_mom_ = set_default_option(options_mom_,'dr_logarithmic_reduction',false); % use logarithmic reduction algorithm to solve the polynomial equation for retrieving the coefficients associated to the endogenous variables in the decision rule
options_mom_ = set_default_option(options_mom_,'dr_logarithmic_reduction_maxiter',100); % maximum number of iterations used in the logarithmic reduction algorithm
options_mom_ = set_default_option(options_mom_,'dr_logarithmic_reduction_tol',1e-12); % convergence criterion used in the cycle reduction algorithm
options_mom_ = set_default_option(options_mom_,'dr_logarithmic_reduction_tol',1e-12); % convergence criterion used in the logarithmic reduction algorithm
options_mom_ = set_default_option(options_mom_,'qz_criterium',1-1e-6); % value used to split stable from unstable eigenvalues in reordering the Generalized Schur decomposition used for solving first order problems
% if there are no unit roots one can use 1.0 (or slightly below) which we set as default; if they are possible, you may have have multiple unit roots and the accuracy decreases when computing the eigenvalues in lyapunov_symm
% Note that unit roots are only possible at first-order, at higher order we set it to 1 in pruned_state_space_system and focus only on stationary observables.
......@@ -453,8 +454,8 @@ if strcmp(mom_method,'IRF_MATCHING') && do_bayesian_estimation
warning('method_of_moments: You specified mh_tune_jscale, but the maximum number of iterations is smaller than the step size. No update will take place.')
end
if options_mom_.load_results_after_load_mh
if ~exist([options_mom_.dirname filesep 'method_of_moments' filesep M_.fname '_mom_results.mat'],'file')
fprintf('\nYou specified the ''load_results_after_load_mh'' option, but no ''%s_mom_results.mat'' file\n',M_.fname);
if ~exist([options_mom_.dirname filesep 'method_of_moments' filesep fname '_mom_results.mat'],'file')
fprintf('\nYou specified the ''load_results_after_load_mh'' option, but no ''%s_mom_results.mat'' file\n',fname);
fprintf('was found in the folder %s%smethod_of_moments.\n',options_mom_.dirname,filesep);
fprintf('Results will be recomputed and option ''load_results_after_load_mh'' is reset to false.\n');
options_mom_.load_results_after_load_mh = false;
......
matlab/+mom/graph_comparison_irfs.m
View file @ 9832f012
......@@ -108,8 +108,9 @@ for jexo = unique_shock_entries' % loop over cell with shock names
% Adding a legend at the bottom
axes('Position',[0, 0, 1, 1],'Visible','off');
lgd = legend([plt_data,plt_model],{'Data', 'Model'}, 'Location', 'southeast','NumColumns',2,'FontSize',14);
lgd.Position = [0.37 0.01 lgd.Position(3) lgd.Position(4)];
if ~isoctave
lgd.Position = [0.37 0.01 lgd.Position(3) lgd.Position(4)];
end
dyn_saveas(fig_irf,[graph_directory_name filesep fname '_matched_irf_' jexo{:} int2str(fig)],nodisplay,graph_format);
if TeX && any(strcmp('eps',cellstr(graph_format)))
fprintf(fid_TeX,'\\begin{figure}[H]\n');
......
matlab/+mom/run.m
View file @ 9832f012
......@@ -95,7 +95,6 @@ function [oo_, options_mom_, M_] = run(bayestopt_, options_, oo_, estim_params_,
% o skipline
% o test_for_deep_parameters_calibration
% o transform_prior_to_laplace_prior
% o warning_config
% Copyright © 2020-2023 Dynare Team
%
......@@ -162,7 +161,7 @@ check_varobs_are_endo_and_declared_once(options_.varobs,M_.endo_names);
% The idea is to be independent of options_ and have full control of the
% estimation instead of possibly having to deal with options chosen somewhere
% else in the mod file.
options_mom_ = mom.default_option_mom_values(options_mom_, options_, M_.dname, do_bayesian_estimation);
options_mom_ = mom.default_option_mom_values(options_mom_, options_, M_.dname, M_.fname, do_bayesian_estimation);
% -------------------------------------------------------------------------
......@@ -182,6 +181,8 @@ end
% -------------------------------------------------------------------------
% initializations
% -------------------------------------------------------------------------
%save warning state for restoring later on
orig_warning_state = warning;
% create output directories to store results
M_.dname = options_mom_.dirname;
CheckPath(M_.dname,'.');
......@@ -689,6 +690,7 @@ if do_bayesian_estimation_mcmc
% skip optimizer-based mode-finding and instead compute the mode based on a run of a MCMC
[~,~,posterior_mode,~] = compute_mh_covariance_matrix(bayestopt_,M_.fname,M_.dname,'method_of_moments');
oo_.mom = fill_mh_mode(posterior_mode',NaN(length(posterior_mode),1),M_,options_mom_,estim_params_,bayestopt_,oo_.mom,'posterior');
warning(orig_warning_state);
return
else
% get stored results if required
......@@ -812,7 +814,7 @@ fprintf('\n==== Method of Moments Estimation (%s) Completed ====\n\n',options_mo
% -------------------------------------------------------------------------
% clean up
% -------------------------------------------------------------------------
warning_config; %reset warning state
warning(orig_warning_state); %reset warning state
if isoctave && isfield(options_mom_, 'prior_restrictions') && ...
isfield(options_mom_.prior_restrictions, 'routine')
% Octave crashes if it tries to save function handles (to the _results.mat file)
......
matlab/+occbin/DSGE_smoother.m
View file @ 9832f012
......@@ -124,6 +124,8 @@ occbin_options.opts_regime.binding_indicator = options_.occbin.smoother.init_bin
occbin_options.opts_regime.regime_history=options_.occbin.smoother.init_regime_history;
error_indicator=false;
options_.noprint = true;
try
%blanket try-catch should be replaced be proper error handling, see https://git.dynare.org/Dynare/dynare/-/merge_requests/2226#note_20318
[alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T0,R0,P,PK,decomp,Trend,state_uncertainty,oo_,bayestopt_,alphahat0,state_uncertainty0] = DsgeSmoother(xparam1,gend,Y,data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_,occbin_options);% T1=TT;
......@@ -178,7 +180,6 @@ opts_simul.endo_init = alphahat0(oo_.dr.inv_order_var,1);
opts_simul.init_regime=regime_history; % use realtime regime for guess, to avoid multiple solution issues!
opts_simul.periods = size(opts_simul.SHOCKS,1);
options_.occbin.simul=opts_simul;
options_.noprint = true;
[~, out, ss] = occbin.solver(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
if out.error_flag
fprintf('Occbin smoother:: simulation within smoother did not converge.\n')
......@@ -479,7 +480,7 @@ if (~is_changed || occbin_smoother_debug) && nargin==12
fprintf(fidTeX,'\\begin{figure}[H]\n');
fprintf(fidTeX,'\\centering \n');
fprintf(fidTeX,'\\includegraphics[width=%2.2f\\textwidth]{%s_smoothedshocks_occbin%s}\n',options_.figures.textwidth*min(j1/3,1),[GraphDirectoryName '/' M_.fname],int2str(ifig)); % don't use filesep as it will create issues with LaTeX on Windows
fprintf(fidTeX,'\\caption{Check plots.}');
fprintf(fidTeX,'\\caption{OccBin smoothed shocks.}');
fprintf(fidTeX,'\\label{Fig:smoothedshocks_occbin:%s}\n',int2str(ifig));
fprintf(fidTeX,'\\end{figure}\n');
fprintf(fidTeX,' \n');
......@@ -488,7 +489,7 @@ if (~is_changed || occbin_smoother_debug) && nargin==12
end
end
if mod(j1,9)~=0 && j==M_.exo_nbr
if mod(j1,9)~=0 && j==M_.exo_nbr
annotation('textbox', [0.1,0,0.35,0.05],'String', 'Linear','Color','Blue','horizontalalignment','center','interpreter','none');
annotation('textbox', [0.55,0,0.35,0.05],'String', 'Piecewise','Color','Red','horizontalalignment','center','interpreter','none');
dyn_saveas(hh_fig,[GraphDirectoryName filesep M_.fname,'_smoothedshocks_occbin',int2str(ifig)],options_.nodisplay,options_.graph_format);
......@@ -497,7 +498,7 @@ if (~is_changed || occbin_smoother_debug) && nargin==12
fprintf(fidTeX,'\\begin{figure}[H]\n');
fprintf(fidTeX,'\\centering \n');
fprintf(fidTeX,'\\includegraphics[width=%2.2f\\textwidth]{%s_smoothedshocks_occbin%s}\n',options_.figures.textwidth*min(j1/3,1),[GraphDirectoryName '/' M_.fname],int2str(ifig)); % don't use filesep as it will create issues with LaTeX on Windows
fprintf(fidTeX,'\\caption{Check plots.}');
fprintf(fidTeX,'\\caption{OccBin smoothed shocks.}');
fprintf(fidTeX,'\\label{Fig:smoothedshocks_occbin:%s}\n',int2str(ifig));
fprintf(fidTeX,'\\end{figure}\n');
fprintf(fidTeX,' \n');
......@@ -505,6 +506,6 @@ if (~is_changed || occbin_smoother_debug) && nargin==12
end
if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
fclose(fidTeX);
end
end
end
end
matlab/+occbin/forecast.m
View file @ 9832f012
function [oo_, error_flag] = forecast(options_,M_,oo_,forecast) %,hist_period)
%function oo_ = forecast(options_,M_,oo_,forecast)
% forecast
function [forecast, error_flag] = forecast(options_,M_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state,forecast_horizon)
% [forecast, error_flag] = forecast(options_,M_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state,forecast_horizon)
% Occbin forecasts
%
% INPUTS
% - options_ [structure] Matlab's structure describing the current options
% - M_ [structure] Matlab's structure describing the model
% - dr_in [structure] model information structure
% - endo_steady_state [double] steady state value for endogenous variables
% - exo_steady_state [double] steady state value for exogenous variables
% - exo_det_steady_state [double] steady state value for exogenous deterministic variables
% - forecast_horizon [integer] forecast horizon
%
% OUTPUTS
% - forecast [structure] forecast results
% - error_flag [integer] error code
%
% SPECIAL REQUIREMENTS
% none.
% Copyright © 2022-2024 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
opts = options_.occbin.forecast;
options_.occbin.simul.maxit = opts.maxit;
options_.occbin.simul.check_ahead_periods = opts.check_ahead_periods;
options_.occbin.simul.periods = forecast;
SHOCKS0 = opts.SHOCKS0;
if ~isempty(SHOCKS0)
if iscell(SHOCKS0)
for j=1:length(SHOCKS0)
sname = SHOCKS0{j}{1};
inds(j)=strmatch(sname,M_.exo_names);
SHOCKS1(j,:)=SHOCKS0{j}{2};
options_.occbin.simul.periods = forecast_horizon;
shocks_input = opts.SHOCKS0;
if ~isempty(shocks_input)
n_shocks=size(shocks_input,1);
if iscell(shocks_input)
inds=NaN(n_shocks,1);
periods=length(shocks_input{1}{2});
shock_mat=NaN(n_shocks,periods);
for j=1:n_shocks
exo_pos=strmatch(shocks_input{j}{1},M_.exo_names,'exact');
if isempty(exo_pos)
error('occbin.forecast: unknown exogenous shock %s',shocks_input{j}{1})
else
inds(j)=exo_pos;
end
if length(shocks_input{j}{2})~=periods
error('occbin.forecast: unknown exogenous shock %s',shocks_input{j}{1})
else
shock_mat(j,:)=shocks_input{j}{2};
end
end
elseif isreal(SHOCKS0)
SHOCKS1=SHOCKS0;
inds = 1:M_.exo_nbr;
elseif isreal(shocks_input)
shock_mat=shocks_input;
inds = (1:M_.exo_nbr)';
end
end
options_.occbin.simul.endo_init = M_.endo_histval(:,1)-endo_steady_state; %initial condition
options_.occbin.simul.init_regime = opts.frcst_regimes;
options_.occbin.simul.init_binding_indicator = [];
shocks_base = zeros(forecast_horizon,M_.exo_nbr);
if ~isempty(shocks_input)
for j=1:n_shocks
shocks_base(:,inds(j))=shock_mat(j,:);
end
end
if opts.replic
options_.noprint=true;
h = dyn_waitbar(0,'Please wait occbin forecast replic ...');
ishock = find(sqrt(diag((M_.Sigma_e))));
options_.occbin.simul.exo_pos=ishock;
effective_exo_nbr= length(ishock);
effective_exo_names = M_.exo_names(ishock);
effective_Sigma_e = M_.Sigma_e(ishock,ishock);
effective_Sigma_e = M_.Sigma_e(ishock,ishock); % does not take heteroskedastic shocks into account
[U,S] = svd(effective_Sigma_e);
% draw random shocks
if opts.qmc
opts.replic =2^(round(log2(opts.replic+1)))-1;
SHOCKS_ant = qmc_sequence(forecast*effective_exo_nbr, int64(1), 1, opts.replic)';
SHOCKS_add = qmc_sequence(forecast_horizon*effective_exo_nbr, int64(1), 1, opts.replic);
else
SHOCKS_ant = randn(forecast*effective_exo_nbr,opts.replic)';
SHOCKS_add = randn(forecast_horizon*effective_exo_nbr,opts.replic);
end
zlin0=zeros(forecast,M_.endo_nbr,opts.replic);
zpiece0=zeros(forecast,M_.endo_nbr,opts.replic);
SHOCKS_add=reshape(SHOCKS_add,effective_exo_nbr,forecast_horizon,opts.replic);
z.linear=NaN(forecast_horizon,M_.endo_nbr,opts.replic);
z.piecewise=NaN(forecast_horizon,M_.endo_nbr,opts.replic);
error_flag=true(opts.replic,1);
simul_SHOCKS=NaN(forecast_horizon,M_.exo_nbr,opts.replic);
for iter=1:opts.replic
if ~isempty(SHOCKS0)
for j=1:length(SHOCKS0)
SHOCKS(:,inds(j))=SHOCKS1(j,:);
end
end
error_flagx=1;
% while error_flagx,
% SHOCKS=transpose(sqrt(diag(diag(effective_Sigma_e)))*(reshape(SHOCKS_ant(iter,:),forecast,effective_exo_nbr))');
SHOCKS=transpose(U*sqrt(S)*(reshape(SHOCKS_ant(iter,:),forecast,effective_exo_nbr))');
% SHOCKS=transpose(U*sqrt(S)*randn(forecast,M_.exo_nbr)'); %realized shocks
options_.occbin.simul.endo_init = M_.endo_histval(:,1)-oo_.steady_state;
options_.occbin.simul.init_regime = opts.frcst_regimes;
options_.occbin.simul.init_binding_indicator = [];
options_.occbin.simul.exo_pos=ishock;
options_.occbin.simul.SHOCKS = SHOCKS;
options_.occbin.simul.SHOCKS = shocks_base(:,ishock)+transpose(U*sqrt(S)*SHOCKS_add(:,:,iter));
options_.occbin.simul.waitbar=0;
[~, out] = occbin.solver(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
zlin0(:,:,iter)=out.linear;
zpiece0(:,:,iter)=out.piecewise;
ys=out.ys;
frcst_regime_history(iter,:)=out.regime_history;
[~, out] = occbin.solver(M_,options_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state);
error_flag(iter)=out.error_flag;
error_flagx = error_flag(iter);
% end
simul_SHOCKS(:,:,iter) = SHOCKS;
if error_flag(iter) && debug_flag
% display('no solution')
% keyboard;
save no_solution SHOCKS zlin0 zpiece0 iter frcst_regime_history
if ~error_flag(iter)
z.linear(:,:,iter)=out.linear;
z.piecewise(:,:,iter)=out.piecewise;
frcst_regime_history(iter,:)=out.regime_history;
error_flag(iter)=out.error_flag;
simul_SHOCKS(:,:,iter) = shocks_base;
else
if options_.debug
save('Occbin_forecast_debug','simul_SHOCKS','z','iter','frcst_regime_history','error_flag','out','shocks_base')
end
end
dyn_waitbar(iter/opts.replic,h,['OccBin MC forecast replic ',int2str(iter),'/',int2str(opts.replic)])
end
dyn_waitbar_close(h);
save temp zlin0 zpiece0 simul_SHOCKS error_flag
if options_.debug
save('Occbin_forecast_debug','simul_SHOCKS','z','iter','frcst_regime_history','error_flag')
end
inx=find(error_flag==0);
zlin0=zlin0(:,:,inx);
zpiece0=zpiece0(:,:,inx);
zlin = mean(zlin0,3);
zpiece = mean(zpiece0,3);
zpiecemin = min(zpiece0,[],3);
zpiecemax = max(zpiece0,[],3);
zlinmin = min(zlin0,[],3);
zlinmax = max(zlin0,[],3);
for i=1:M_.endo_nbr
for j=1:forecast
[post_mean(j,1), post_median(j,1), post_var(j,1), hpd_interval(j,:), post_deciles(j,:)] = posterior_moments(squeeze(zlin0(j,i,:)),options_.forecasts.conf_sig);
end
oo_.occbin.linear_forecast.Mean.(M_.endo_names{i})=post_mean;
oo_.occbin.linear_forecast.Median.(M_.endo_names{i})=post_median;
oo_.occbin.linear_forecast.Var.(M_.endo_names{i})=post_var;
oo_.occbin.linear_forecast.HPDinf.(M_.endo_names{i})=hpd_interval(:,1);
oo_.occbin.linear_forecast.HPDsup.(M_.endo_names{i})=hpd_interval(:,2);
oo_.occbin.linear_forecast.Deciles.(M_.endo_names{i})=post_deciles;
oo_.occbin.linear_forecast.Min.(M_.endo_names{i})=zlinmin(:,i);
oo_.occbin.linear_forecast.Max.(M_.endo_names{i})=zlinmax(:,i);
for j=1:forecast
[post_mean(j,1), post_median(j,1), post_var(j,1), hpd_interval(j,:), post_deciles(j,:)] = posterior_moments(squeeze(zpiece0(j,i,:)),options_.forecasts.conf_sig);
end
oo_.occbin.forecast.Mean.(M_.endo_names{i})=post_mean;
oo_.occbin.forecast.Median.(M_.endo_names{i})=post_median;
oo_.occbin.forecast.Var.(M_.endo_names{i})=post_var;
oo_.occbin.forecast.HPDinf.(M_.endo_names{i})=hpd_interval(:,1);
oo_.occbin.forecast.HPDsup.(M_.endo_names{i})=hpd_interval(:,2);
oo_.occbin.forecast.Deciles.(M_.endo_names{i})=post_deciles;
oo_.occbin.forecast.Min.(M_.endo_names{i})=zpiecemin(:,i);
oo_.occbin.forecast.Max.(M_.endo_names{i})=zpiecemax(:,i);
% eval([M_.endo_names{i},'_ss=zdatass(i);']);
if length(inx)<opts.replic
fprintf('\noccbin.forecast: forecast simulation was only successful in %u of %u simulations.\n',length(inx),opts.replic);
end
else
SHOCKS = zeros(forecast,M_.exo_nbr);
if ~isempty(SHOCKS0)
for j=1:length(SHOCKS0)
SHOCKS(:,inds(j))=SHOCKS1(j,:);
z.linear=z.linear(:,:,inx);
z.piecewise=z.piecewise(:,:,inx);
z.min.piecewise = min(z.piecewise,[],3);
z.max.piecewise = max(z.piecewise,[],3);
z.min.linear = min(z.linear,[],3);
z.max.linear = max(z.linear,[],3);
field_names={'linear','piecewise'};
post_mean=NaN(forecast_horizon,1);
post_median=NaN(forecast_horizon,1);
post_var=NaN(forecast_horizon,1);
hpd_interval=NaN(forecast_horizon,2);
post_deciles=NaN(forecast_horizon,9);
for field_iter=1:2
for i=1:M_.endo_nbr
for j=1:forecast_horizon
[post_mean(j,1), post_median(j,1), post_var(j,1), hpd_interval(j,:), post_deciles(j,:)] = posterior_moments(squeeze(z.(field_names{field_iter})(j,i,:)),options_.forecasts.conf_sig);
end
forecast.(field_names{field_iter}).Mean.(M_.endo_names{i})=post_mean;
forecast.(field_names{field_iter}).Median.(M_.endo_names{i})=post_median;
forecast.(field_names{field_iter}).Var.(M_.endo_names{i})=post_var;
forecast.(field_names{field_iter}).HPDinf.(M_.endo_names{i})=hpd_interval(:,1);
forecast.(field_names{field_iter}).HPDsup.(M_.endo_names{i})=hpd_interval(:,2);
forecast.(field_names{field_iter}).Deciles.(M_.endo_names{i})=post_deciles;
forecast.(field_names{field_iter}).Min.(M_.endo_names{i})=z.min.(field_names{field_iter})(:,i);
forecast.(field_names{field_iter}).Max.(M_.endo_names{i})=z.max.(field_names{field_iter})(:,i);
end
end
options_.occbin.simul.endo_init = M_.endo_histval(:,1)-oo_.steady_state;
options_.occbin.simul.init_regime = opts.frcst_regimes;
options_.occbin.simul.init_violvecbool = [];
else
options_.occbin.simul.irfshock = M_.exo_names;
options_.occbin.simul.SHOCKS = SHOCKS;
[~, out] = occbin.solver(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
options_.occbin.simul.SHOCKS = shocks_base;
[~, out] = occbin.solver(M_,options_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state);
error_flag=out.error_flag;
if out.error_flag
fprintf('occbin.forecast: forecast simulation failed.')
return;
end
zlin=out.linear;
zpiece=out.piecewise;
frcst_regime_history=out.regime_history;
error_flag=out.error_flag;
for i=1:M_.endo_nbr
oo_.occbin.linear_forecast.Mean.(M_.endo_names{i})= zlin(:,i);
oo_.occbin.forecast.Mean.(M_.endo_names{i})= zpiece(:,i);
oo_.occbin.forecast.HPDinf.(M_.endo_names{i})= nan;
oo_.occbin.forecast.HPDsup.(M_.endo_names{i})= nan;
forecast.linear.Mean.(M_.endo_names{i})= out.linear(:,i);
forecast.piecewise.Mean.(M_.endo_names{i})= out.piecewise(:,i);
end
end
oo_.occbin.forecast.regimes=frcst_regime_history;
forecast.regimes=frcst_regime_history;
\ No newline at end of file
matlab/+occbin/irf.m
View file @ 9832f012
function [oo_] = irf(M_,oo_,options_)
function irfs = irf(M_,oo_,options_)
% irfs = irf(M_,oo_,options_)
% Calls a minimizer
%
% INPUTS
% - M_ [structure] Matlab's structure describing the model
% - oo_ [structure] Matlab's structure containing the results
% - options_ [structure] Matlab's structure describing the current options
%
% OUTPUTS
% - irfs [structure] IRF results
%
% SPECIAL REQUIREMENTS
% none.
%
%
% Copyright © 2022-2023 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
shocknames = options_.occbin.irf.exo_names;
shocksigns = options_.occbin.irf.shocksigns;
shocksigns = options_.occbin.irf.shocksigns; %'pos','neg'
shocksize = options_.occbin.irf.shocksize;
t0 = options_.occbin.irf.t0;
t_0 = options_.occbin.irf.t0;
%% set simulation options based on IRF options
options_.occbin.simul.init_regime = options_.occbin.irf.init_regime;
options_.occbin.simul.check_ahead_periods = options_.occbin.irf.check_ahead_periods;
options_.occbin.simul.maxit = options_.occbin.irf.maxit;
options_.occbin.simul.periods = options_.irf;
% Run inital conditions + other shocks
if t0 == 0
shocks0= zeros(options_.occbin.simul.periods+1,M_.exo_nbr);
%% Run initial conditions + other shocks
if t_0 == 0
shocks_base = zeros(options_.occbin.simul.periods+1,M_.exo_nbr);
options_.occbin.simul.endo_init = [];
else
% girf conditional to smoothed states in t0 and shocks in t0+1
shocks0= [oo_.occbin.smoother.etahat(:,t0+1)'; zeros(options_.occbin.simul.periods,M_.exo_nbr)];
options_.occbin.simul.SHOCKS=shocks0;
options_.occbin.simul.endo_init = oo_.occbin.smoother.alphahat(oo_.dr.inv_order_var,t0);
if ~isfield(oo_.occbin,'smoother')
error('occbin.irfs: smoother must be run before requesting GIRFs based on smoothed results')
end
% GIRF conditional on smoothed states in t_0 and shocks in t_0+1
shocks_base= [oo_.occbin.smoother.etahat(:,t_0+1)'; zeros(options_.occbin.simul.periods,M_.exo_nbr)];
options_.occbin.simul.SHOCKS=shocks_base;
options_.occbin.simul.endo_init = oo_.occbin.smoother.alphahat(oo_.dr.inv_order_var,t_0);
end
options_.occbin.simul.SHOCKS=shocks_base;
[~, out_base] = occbin.solver(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
if out_base.error_flag
error('occbin.irfs: could not compute the solution')
end
options_.occbin.simul.SHOCKS=shocks0;
[~, out0] = occbin.solver(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
zlin0 = out0.linear;
zpiece0 = out0.piecewise;
% Select shocks of interest
jexo_all =zeros(size(shocknames,1),1);
irfs.linear = struct();
irfs.piecewise = struct();
% Get indices of shocks of interest
exo_index =zeros(size(shocknames,1),1);
for i=1:length(shocknames)
jexo_all(i) = strmatch(shocknames{i},M_.exo_names,'exact');
exo_index(i) = strmatch(shocknames{i},M_.exo_names,'exact');
end
oo_.occbin.linear_irfs = struct();
oo_.occbin.irfs = struct();
% cs=get_lower_cholesky_covariance(M_.Sigma_e,options_.add_tiny_number_to_cholesky);
% irf_shocks_indx = getIrfShocksIndx(M_, options_);
% Set shock size
if isempty(shocksize)
% if isfield(oo_.posterior_mode.shocks_std,M_.exo_names{jexo})
shocksize = sqrt(diag(M_.Sigma_e(jexo_all,jexo_all))); %oo_.posterior_mode.shocks_std.(M_.exo_names{jexo});
shocksize = sqrt(diag(M_.Sigma_e(exo_index,exo_index)));
if any(shocksize < 1.e-9)
shocksize(shocksize < 1.e-9) = 0.01;
end
end
if numel(shocksize)==1
shocksize=repmat(shocksize,[length(shocknames),1]);
end
% Run IRFs
for counter = 1:length(jexo_all)
jexo = jexo_all(counter);
if ~options_.noprint
fprintf('Producing GIRFs for shock %s. Simulation %d out of %d. \n',M_.exo_names{jexo},counter,size(jexo_all,1));
end
if ismember('pos',shocksigns)
% (+) shock
shocks1=shocks0;
shocks1(1,jexo)=shocks0(1,jexo)+shocksize(counter);
if t0 == 0
options_.occbin.simul.SHOCKS=shocks1;
options_.occbin.simul.endo_init = [];
[~, out_pos] = occbin.solver(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
else
options_.occbin.simul.SHOCKS=shocks1;
options_.occbin.simul.endo_init = oo_.occbin.smoother.alphahat(oo_.dr.inv_order_var,t0);
[~, out_pos] = occbin.solver(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
for sign_iter=1:length(shocksigns)
for IRF_counter = 1:length(exo_index)
jexo = exo_index(IRF_counter);
if ~options_.noprint && options_.debug
fprintf('occbin.irf: Producing GIRFs for shock %s. Simulation %d out of %d. \n',M_.exo_names{jexo},IRF_counter,size(exo_index,1));
end
if out_pos.error_flag
warning('Occbin error.')
return
shocks1=shocks_base;
if ismember('pos',shocksigns{sign_iter})
shocks1(1,jexo)=shocks_base(1,jexo)+shocksize(IRF_counter);
elseif ismember('neg',shocksigns{sign_iter})
shocks1(1,jexo)=shocks_base(1,jexo)-shocksize(IRF_counter);
end
zlin_pos = out_pos.linear;
zpiece_pos = out_pos.piecewise;
% Substract inital conditions + other shocks
zlin_pos_diff = zlin_pos-zlin0;
zpiece_pos_diff = zpiece_pos-zpiece0;
end
if ismember('neg',shocksigns)
% (-) shock
shocks_1=shocks0;
shocks_1(1,jexo)=shocks0(1,jexo)-shocksize(counter);
if t0 == 0
options_.occbin.simul.SHOCKS=shocks_1;
options_.occbin.simul.SHOCKS=shocks1;
if t_0 == 0
options_.occbin.simul.endo_init = [];
[~, out_neg] = occbin.solver(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
else
options_.occbin.simul.SHOCKS=shocks_1;
options_.occbin.simul.endo_init = oo_.occbin.smoother.alphahat(oo_.dr.inv_order_var,t0);
[~, out_neg] = occbin.solver(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
options_.occbin.simul.endo_init = oo_.occbin.smoother.alphahat(oo_.dr.inv_order_var,t_0);
end
if out_neg.error_flag
warning('Occbin error.')
return
[~, out_sim] = occbin.solver(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
if out_sim.error_flag
warning('occbin.irfs: simulation failed')
skip
end
zlin_neg = out_neg.linear;
zpiece_neg = out_neg.piecewise;
zlin_neg_diff = zlin_neg-zlin0;
zpiece_neg_diff = zpiece_neg-zpiece0;
end
% Save
if ~isfield(oo_.occbin,'linear_irfs')
oo_.occbin.linear_irfs = struct();
end
if ~isfield(oo_.occbin,'irfs')
oo_.occbin.irfs = struct();
end
% Substract inital conditions + other shocks
zdiff.linear.(shocksigns{sign_iter}) = out_sim.linear-out_base.linear;
zdiff.piecewise.(shocksigns{sign_iter}) = out_sim.piecewise-out_base.piecewise;
for jendo=1:M_.endo_nbr
% oo_.occbin.irfs.([M_.endo_names{jendo} '_' M_.exo_names{jexo} '1']) = zpiece_pos (:,jendo);
% oo_.occbin.irfs.([M_.endo_names{jendo} '_' M_.exo_names{jexo} '_1']) = zpiece_neg (:,jendo);
% oo_.occbin.linear_irfs.([M_.endo_names{jendo} '_' M_.exo_names{jexo} '1']) = zlin_pos (:,jendo);
% oo_.occbin.linear_irfs.([M_.endo_names{jendo} '_' M_.exo_names{jexo} '_1']) = zlin_neg(:,jendo);
if ismember('pos',shocksigns)
oo_.occbin.irfs.([M_.endo_names{jendo} '_' M_.exo_names{jexo} '_pos']) = zpiece_pos_diff (:,jendo);
oo_.occbin.linear_irfs.([M_.endo_names{jendo} '_' M_.exo_names{jexo} '_pos']) = zlin_pos_diff (:,jendo);
for j_endo=1:M_.endo_nbr
if ismember('pos',shocksigns)
irfs.piecewise.([M_.endo_names{j_endo} '_' M_.exo_names{jexo} '_' shocksigns{sign_iter}]) = zdiff.piecewise.(shocksigns{sign_iter})(:,j_endo);
irfs.linear.([M_.endo_names{j_endo} '_' M_.exo_names{jexo} '_' shocksigns{sign_iter}]) = zdiff.linear.(shocksigns{sign_iter})(:,j_endo);
end
end
if ismember('neg',shocksigns)
oo_.occbin.irfs.([M_.endo_names{jendo} '_' M_.exo_names{jexo} '_neg']) = zpiece_neg_diff (:,jendo);
oo_.occbin.linear_irfs.([M_.endo_names{jendo} '_' M_.exo_names{jexo} '_neg']) = zlin_neg_diff (:,jendo);
end
% %
% oo_.occbin.irfs0.([M_.endo_names{jendo} '_' M_.exo_names{jexo} '1']) = zpiece0(:,jendo);
% oo_.occbin.linear_irfs0.([M_.endo_names{jendo} '_' M_.exo_names{jexo} '1']) = zlin0(:,jendo);
% oo_.occbin.irfs0.([M_.endo_names{jendo} '_' M_.exo_names{jexo} '_1']) = zpiece0(:,jendo);
% oo_.occbin.linear_irfs0.([M_.endo_names{jendo} '_' M_.exo_names{jexo} '_1']) = zlin0(:,jendo);
end
end
end
end
\ No newline at end of file
matlab/+occbin/kalman_update_algo_1.m
View file @ 9832f012
......@@ -62,7 +62,7 @@ function [a, a1, P, P1, v, T, R, C, regimes_, error_flag, M_, lik, etahat, alpha
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
warning off
orig_warning_state = warning;
options_.noprint = true;
R=NaN(size(RR));
......@@ -148,6 +148,7 @@ else
end
if error_flag
etahat=NaN(size(QQQ,1),1);
warning(orig_warning_state);
return;
end
......@@ -183,6 +184,7 @@ if out.error_flag
error_flag = out.error_flag;
etahat=etahat(:,2);
lik=inf;
warning(orig_warning_state);
return;
end
......@@ -211,7 +213,12 @@ if any(myregime) || ~isequal(regimes_(1),regimes0(1))
if M_.occbin.constraint_nbr==1
regime_end(niter) = regimes_(1).regimestart(end);
end
[a, a1, P, P1, v, alphahat, etahat, lik, V] = occbin_kalman_update0(a,a1,P,P1,data_index,Z,v,Y,H,QQQ,TT,RR,CC,iF,L,mm, options_.rescale_prediction_error_covariance, options_.occbin.likelihood.IF_likelihood, options_.occbin.filter.state_covariance);
[a, a1, P, P1, v, alphahat, etahat, lik, V, error_flag] = occbin_kalman_update0(a,a1,P,P1,data_index,Z,v,Y,H,QQQ,TT,RR,CC,iF,L,mm, options_.rescale_prediction_error_covariance, options_.occbin.likelihood.IF_likelihood, options_.occbin.filter.state_covariance);
if error_flag
etahat=NaN(size(QQQ,1),1);
warning(orig_warning_state);
return;
end
etahat_hist(niter) = {etahat};
lik_hist(niter) = lik;
opts_simul.SHOCKS(1,:) = etahat(:,2)';
......@@ -239,6 +246,7 @@ if any(myregime) || ~isequal(regimes_(1),regimes0(1))
error_flag = out.error_flag;
etahat=etahat(:,2);
lik=inf;
warning(orig_warning_state);
return;
end
regimes0=regimes_;
......@@ -268,18 +276,25 @@ if any(myregime) || ~isequal(regimes_(1),regimes0(1))
error_flag = out.error_flag;
etahat=etahat(:,2);
lik=inf;
warning(orig_warning_state);
return;
else
regimes_ = out.regime_history;
TT(:,:,2)=ss.T(my_order_var,my_order_var,1);
RR(:,:,2)=ss.R(my_order_var,:,1);
CC(:,2)=ss.C(my_order_var,1);
[a, a1, P, P1, v, alphahat, etahat, lik, V] = occbin_kalman_update0(a,a1,P,P1,data_index,Z,v,Y,H,QQQ,TT,RR,CC,iF,L,mm, options_.rescale_prediction_error_covariance, options_.occbin.likelihood.IF_likelihood, options_.occbin.filter.state_covariance);
[a, a1, P, P1, v, alphahat, etahat, lik, V, error_flag] = occbin_kalman_update0(a,a1,P,P1,data_index,Z,v,Y,H,QQQ,TT,RR,CC,iF,L,mm, options_.rescale_prediction_error_covariance, options_.occbin.likelihood.IF_likelihood, options_.occbin.filter.state_covariance);
if error_flag
etahat=NaN(size(QQQ,1),1);
warning(orig_warning_state);
return;
end
end
else
error_flag = 330;
etahat=etahat(:,2);
lik=inf;
warning(orig_warning_state);
return;
end
end
......@@ -289,7 +304,7 @@ end
error_flag = out.error_flag;
if ~error_flag && niter>options_.occbin.likelihood.max_number_of_iterations && ~isequal(regimes_(1),regimes0(1))
error_flag = 1;
error_flag = 331;
end
if ~error_flag
......@@ -306,7 +321,8 @@ if nargout>=13
etahat=etahat(:,2);
end
warning_config;
%reset warning state
warning(orig_warning_state);
end
function [a, a1, P, P1, v, alphahat, etahat, lik, V, error_flag] = occbin_kalman_update0(a,a1,P,P1,data_index,Z,v,Y,H,QQQ,TT,RR,CC,iF,L,mm, rescale_prediction_error_covariance, IF_likelihood, state_uncertainty_flag)
......@@ -322,7 +338,7 @@ else
V=[];
end
warning off
orig_warning_state = warning;
if nargin<18
IF_likelihood=0;
end
......@@ -351,7 +367,8 @@ else
F = ZZ*P(:,:,t)*ZZ' + H(di,di);
sig=sqrt(diag(F));
if any(any(isnan(F)))
error_flag=1;
error_flag=325;
warning(orig_warning_state);
return;
end
if rank(F)<size(F,1)
......@@ -434,5 +451,6 @@ while t > 1
end
end
warning_config;
%reset warning state
warning(orig_warning_state);
end