From c90368d48cac1a6cb9dc42fb1cb32714f32846fc Mon Sep 17 00:00:00 2001 From: Johannes Pfeifer Date: Wed, 23 Nov 2022 15:17:51 +0100 Subject: [PATCH 1/6] NKM.mod: clean up file --- tests/occbin/filter/NKM.mod | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/tests/occbin/filter/NKM.mod b/tests/occbin/filter/NKM.mod index 8f902e1d4..8c04a966e 100644 --- a/tests/occbin/filter/NKM.mod +++ b/tests/occbin/filter/NKM.mod @@ -11,10 +11,6 @@ @#define small_model = 0 @#endif -// if ~exist('run_ivf','var') -run_ivf=0; -// end - // ----------------- Defintions -----------------------------------------// var c //1 Consumption @@ -278,8 +274,6 @@ check; varobs yg inom pi; estimated_params; - // PARAM NAME, INITVAL, LB, UB, PRIOR_SHAPE, PRIOR_P1, PRIOR_P2, PRIOR_P3, PRIOR_P4, JSCALE - // PRIOR_SHAPE: BETA_PDF, GAMMA_PDF, NORMAL_PDF, INV_GAMMA_PDF varphip,,0,inf,NORMAL_PDF,100,25; phipi,,,,NORMAL_PDF,2,0.25; phiy,,0,inf,NORMAL_PDF,0.5,0.25; @@ -291,8 +285,6 @@ varobs yg inom pi; sigi,,,,INV_GAMMA_PDF,0.002,0.002; end; - -// dataloading_jme_beta(1,'sims.txt',30); load('dataobsfile','inom') // check if inom is at lb and remove data + associated shock verbatim; @@ -345,4 +337,12 @@ varobs yg inom pi; subplot(223) plot([oo0.SmoothedShocks.epss oo_.SmoothedShocks.epss]), title('epss') legend('PKF','IF') + figure, + subplot(221) + plot([oo0.SmoothedVariables.inom oo_.SmoothedVariables.inom]), title('inom') + subplot(222) + plot([oo0.SmoothedVariables.yg oo_.SmoothedVariables.yg]), title('yg') + subplot(223) + plot([oo0.SmoothedVariables.pi oo_.SmoothedVariables.pi]), title('pi') + legend('PKF','IF') occbin_write_regimes(smoother); From eeecccd29b6cedaf843328d58486826b0d7b49c6 Mon Sep 17 00:00:00 2001 From: Johannes Pfeifer Date: Wed, 23 Nov 2022 15:18:44 +0100 Subject: [PATCH 2/6] IVF_core.m: fix header --- matlab/+occbin/IVF_core.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/matlab/+occbin/IVF_core.m b/matlab/+occbin/IVF_core.m index efc3f6e84..ad943bafb 100644 --- a/matlab/+occbin/IVF_core.m +++ b/matlab/+occbin/IVF_core.m @@ -1,5 +1,5 @@ function [filtered_errs, resids, Emat, stateval, error_code] = IVF_core(M_,oo_,options_,err_index,filtered_errs_init,my_obs_list,obs,init_val) -% function [filtered_errs, resids, Emat, stateval] = IVF_core(M_,oo_,options_,err_index,filtered_errs_init,my_obs_list,obs,init_val) +% function [filtered_errs, resids, Emat, stateval, error_code] = IVF_core(M_,oo_,options_,err_index,filtered_errs_init,my_obs_list,obs,init_val) % Computes thre % % Outputs: From 1e2fb88d326bed478ba9f69e37fc89db628457c6 Mon Sep 17 00:00:00 2001 From: Johannes Pfeifer Date: Wed, 23 Nov 2022 15:37:12 +0100 Subject: [PATCH 3/6] IVF: fix error handling for smoother --- matlab/+occbin/IVF_posterior.m | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/matlab/+occbin/IVF_posterior.m b/matlab/+occbin/IVF_posterior.m index 5959d8ba1..2783e2b1c 100644 --- a/matlab/+occbin/IVF_posterior.m +++ b/matlab/+occbin/IVF_posterior.m @@ -106,7 +106,14 @@ filtered_errs_init = zeros(sample_length,sum(err_index)); if info(1) fval = Inf; exit_flag = 0; + atT=NaN(size(stateval(:,DynareResults.dr.order_var)')); + innov=NaN(Model.exo_nbr,sample_length); return +else + atT = stateval(:,DynareResults.dr.order_var)'; + innov = zeros(Model.exo_nbr,sample_length); + innov(diag(Model.Sigma_e)~=0,:)=filtered_errs'; + updated_variables = atT*nan; end nobs=size(filtered_errs,1); @@ -185,10 +192,7 @@ end % remember that the likelihood has already been multiplied by -1 % hence, posterior is -1 times the log of the prior fval = like+prior; -atT = stateval(:,DynareResults.dr.order_var)'; -innov = zeros(Model.exo_nbr,sample_length); -innov(diag(Model.Sigma_e)~=0,:)=filtered_errs'; -updated_variables = atT*nan; + BayesInfo.mf = BayesInfo.smoother_var_list(BayesInfo.smoother_mf); From 80f6799c38293bdb7d72ac67be071ce3926dd2ee Mon Sep 17 00:00:00 2001 From: Johannes Pfeifer Date: Wed, 23 Nov 2022 16:44:18 +0100 Subject: [PATCH 4/6] IVF: only store smoother results if requested --- matlab/+occbin/IVF_posterior.m | 11 +---------- 1 file changed, 1 insertion(+), 10 deletions(-) diff --git a/matlab/+occbin/IVF_posterior.m b/matlab/+occbin/IVF_posterior.m index 2783e2b1c..db4d51ec0 100644 --- a/matlab/+occbin/IVF_posterior.m +++ b/matlab/+occbin/IVF_posterior.m @@ -113,7 +113,6 @@ else atT = stateval(:,DynareResults.dr.order_var)'; innov = zeros(Model.exo_nbr,sample_length); innov(diag(Model.Sigma_e)~=0,:)=filtered_errs'; - updated_variables = atT*nan; end nobs=size(filtered_errs,1); @@ -191,12 +190,4 @@ end % remember that the likelihood has already been multiplied by -1 % hence, posterior is -1 times the log of the prior -fval = like+prior; - -BayesInfo.mf = BayesInfo.smoother_var_list(BayesInfo.smoother_mf); - - -initDynareOptions=DynareOptions; -DynareOptions.nk=[]; %unset options_.nk and reset it later -[DynareResults]=store_smoother_results(Model,DynareResults,DynareOptions,BayesInfo,dataset_,obs_info,atT,innov,[],updated_variables,DynareResults.dr.ys,zeros(length(DynareOptions.varobs_id),1)); -DynareOptions=initDynareOptions; +fval = like+prior; \ No newline at end of file From 26fbc6c56da1f316653f326c5a004746b3012bb9 Mon Sep 17 00:00:00 2001 From: Johannes Pfeifer Date: Wed, 23 Nov 2022 16:46:40 +0100 Subject: [PATCH 5/6] IVF: improve error handling --- matlab/dynare_estimation_1.m | 42 ++++++++++++++++++++++++------------ matlab/evaluate_smoother.m | 26 ++++++++++++++++------ 2 files changed, 47 insertions(+), 21 deletions(-) diff --git a/matlab/dynare_estimation_1.m b/matlab/dynare_estimation_1.m index 7465276de..9e8d1af5e 100644 --- a/matlab/dynare_estimation_1.m +++ b/matlab/dynare_estimation_1.m @@ -175,13 +175,20 @@ if isequal(options_.mode_compute,0) && isempty(options_.mode_file) && ~options_. if options_.order==1 && ~options_.particle.status if options_.smoother if options_.occbin.smoother.status && options_.occbin.smoother.inversion_filter - [~, ~, ~, ~, ~, ~, ~, ~, ~, ~, oo_, atT, innov] = occbin.IVF_posterior(xparam1,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,prior_bounds(bayestopt_,options_.prior_trunc),oo_); - updated_variables = atT*nan; - measurement_error=[]; - ys = oo_.dr.ys; - trend_coeff = zeros(length(options_.varobs_id),1); - bayestopt_.mf = bayestopt_.smoother_var_list(bayestopt_.smoother_mf); - [oo_]=store_smoother_results(M_,oo_,options_,bayestopt_,dataset_,dataset_info,atT,innov,measurement_error,updated_variables,ys,trend_coeff); + [~, info, ~, ~, ~, ~, ~, ~, ~, ~, oo_, atT, innov] = occbin.IVF_posterior(xparam1,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,prior_bounds(bayestopt_,options_.prior_trunc),oo_); + if ismember(info(1),[303,304,306]) + fprintf('\nIVF: smoother did not succeed. No results will be written to oo_.\n') + else + updated_variables = atT*nan; + measurement_error=[]; + ys = oo_.dr.ys; + trend_coeff = zeros(length(options_.varobs_id),1); + bayestopt_.mf = bayestopt_.smoother_var_list(bayestopt_.smoother_mf); + options_nk=options_.nk; + options_.nk=[]; %unset options_.nk and reset it later + [oo_]=store_smoother_results(M_,oo_,options_,bayestopt_,dataset_,dataset_info,atT,innov,measurement_error,updated_variables,ys,trend_coeff); + options_.nk=options_nk; + end else if options_.occbin.smoother.status [atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,T,R,P,PK,decomp,Trend,state_uncertainty,M_,oo_,bayestopt_] = occbin.DSGE_smoother(xparam1,gend,transpose(data),data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_,dataset_,dataset_info); @@ -593,13 +600,20 @@ if (~((any(bayestopt_.pshape > 0) && options_.mh_replic) || (any(bayestopt_.psha || ~options_.smoother ) && ~options_.partial_information % to be fixed %% ML estimation, or posterior mode without Metropolis-Hastings or Metropolis without Bayesian smoothed variables if options_.occbin.smoother.status && options_.occbin.smoother.inversion_filter - [~, ~, ~, ~, ~, ~, ~, ~, ~, ~, oo_, atT, innov] = occbin.IVF_posterior(xparam1,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,prior_bounds(bayestopt_,options_.prior_trunc),oo_); - updated_variables = atT*nan; - measurement_error=[]; - ys = oo_.dr.ys; - trend_coeff = zeros(length(options_.varobs_id),1); - bayestopt_.mf = bayestopt_.smoother_var_list(bayestopt_.smoother_mf); - [oo_, yf]=store_smoother_results(M_,oo_,options_,bayestopt_,dataset_,dataset_info,atT,innov,measurement_error,updated_variables,ys,trend_coeff); + [~, info, ~, ~, ~, ~, ~, ~, ~, ~, oo_, atT, innov] = occbin.IVF_posterior(xparam1,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,prior_bounds(bayestopt_,options_.prior_trunc),oo_); + if ismember(info(1),[303,304,306]) + fprintf('\nIVF: smoother did not succeed. No results will be written to oo_.\n') + else + updated_variables = atT*nan; + measurement_error=[]; + ys = oo_.dr.ys; + trend_coeff = zeros(length(options_.varobs_id),1); + bayestopt_.mf = bayestopt_.smoother_var_list(bayestopt_.smoother_mf); + options_nk=options_.nk; + options_.nk=[]; %unset options_.nk and reset it later + [oo_, yf]=store_smoother_results(M_,oo_,options_,bayestopt_,dataset_,dataset_info,atT,innov,measurement_error,updated_variables,ys,trend_coeff); + options_.nk=options_nk; + end else if options_.occbin.smoother.status [atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,T,R,P,PK,decomp,Trend,state_uncertainty,M_,oo_,bayestopt_] = occbin.DSGE_smoother(xparam1,dataset_.nobs,transpose(dataset_.data),dataset_info.missing.aindex,dataset_info.missing.state,M_,oo_,options_,bayestopt_,estim_params_,dataset_,dataset_info); diff --git a/matlab/evaluate_smoother.m b/matlab/evaluate_smoother.m index a03cff563..b0fb55501 100644 --- a/matlab/evaluate_smoother.m +++ b/matlab/evaluate_smoother.m @@ -105,12 +105,17 @@ end if options_.occbin.smoother.status if options_.occbin.smoother.inversion_filter - [~, ~, ~, ~, ~, ~, ~, ~, ~, ~, oo_, atT, innov] = occbin.IVF_posterior(parameters,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,prior_bounds(bayestopt_,options_.prior_trunc),oo_); - updated_variables = atT*nan; - measurement_error=[]; - ys = oo_.dr.ys; - trend_coeff = zeros(length(options_.varobs_id),1); - bayestopt_.mf = bayestopt_.smoother_var_list(bayestopt_.smoother_mf); + [~, info, ~, ~, ~, ~, ~, ~, ~, ~, oo_, atT, innov] = occbin.IVF_posterior(parameters,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,prior_bounds(bayestopt_,options_.prior_trunc),oo_); + if ismember(info(1),[303,304,306]) + oo_.occbin.smoother.error_flag=1; + else + oo_.occbin.smoother.error_flag=0; + updated_variables = atT*nan; + measurement_error=[]; + ys = oo_.dr.ys; + trend_coeff = zeros(length(options_.varobs_id),1); + bayestopt_.mf = bayestopt_.smoother_var_list(bayestopt_.smoother_mf); + end else [atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,T,R,P,PK,decomp,Trend,state_uncertainty,M_,oo_,bayestopt_] = ... occbin.DSGE_smoother(parameters,dataset_.nobs,transpose(dataset_.data),dataset_info.missing.aindex,dataset_info.missing.state,M_,oo_,options_,bayestopt_,estim_params_,dataset_,dataset_info); @@ -124,7 +129,14 @@ if ~(options_.occbin.smoother.status && options_.occbin.smoother.inversion_filte [oo_]=store_smoother_results(M_,oo_,options_,bayestopt_,dataset_,dataset_info,atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,P,PK,decomp,Trend,state_uncertainty); end else - [oo_]=store_smoother_results(M_,oo_,options_,bayestopt_,dataset_,dataset_info,atT,innov,measurement_error,updated_variables,ys,trend_coeff); + if ~oo_.occbin.smoother.error_flag + options_nk=options_.nk; + options_.nk=[]; %unset options_.nk and reset it later + [oo_]=store_smoother_results(M_,oo_,options_,bayestopt_,dataset_,dataset_info,atT,innov,measurement_error,updated_variables,ys,trend_coeff); + options_.nk=options_nk; + else + fprintf('\nIVF: smoother did not succeed. No results will be written to oo_.\n') + end end if nargout>4 Smoothed_variables_declaration_order_deviation_form=atT(oo_.dr.inv_order_var(bayestopt_.smoother_var_list),:); From e530f7cbaa83505b8d7acf6844917694857cda89 Mon Sep 17 00:00:00 2001 From: Johannes Pfeifer Date: Wed, 23 Nov 2022 16:48:08 +0100 Subject: [PATCH 6/6] OccBin: add calib_smoother testfile --- tests/Makefile.am | 1 + tests/occbin/filter/NKM_calib_smoother.mod | 334 +++++++++++++++++++++ 2 files changed, 335 insertions(+) create mode 100644 tests/occbin/filter/NKM_calib_smoother.mod diff --git a/tests/Makefile.am b/tests/Makefile.am index 0b2e63789..07d7ec6f8 100644 --- a/tests/Makefile.am +++ b/tests/Makefile.am @@ -5,6 +5,7 @@ MODFILES = \ occbin/model_borrcon/borrcon.mod \ occbin/model_borrcon/borrcon_0_std_shocks.mod \ occbin/filter/NKM.mod \ + occbin/filter/NKM_calib_smoother.mod \ occbin/filter/NKM_0_std_shocks.mod \ optimizers/fs2000_6.mod \ moments/example1_hp_test.mod \ diff --git a/tests/occbin/filter/NKM_calib_smoother.mod b/tests/occbin/filter/NKM_calib_smoother.mod new file mode 100644 index 000000000..2e81f6fa9 --- /dev/null +++ b/tests/occbin/filter/NKM_calib_smoother.mod @@ -0,0 +1,334 @@ +//Tests Occbin estimation with IVF and PKF with 1 constraints + +// this file implements the model in: +// Atkinson, T., A. W. Richter, and N. A. Throckmorton (2019). +// The zero lower bound andestimation accuracy.Journal of Monetary Economics +// original codes provided by Alexander Richter +// adapted for dynare implementation +// ------------------------- Settings -----------------------------------------// + +@#ifndef small_model + @#define small_model = 0 +@#endif + + +// ----------------- Defintions -----------------------------------------// +var + c //1 Consumption + n //2 Labor + y //5 Output + yf //6 Final goods + yg //11 Output growth gap + w //12 Real wage rate + wf //13 Flexible real wage + pigap //15 Inflation rate -> pi(t)/pibar = pigap + inom //16 Nominal interest rate + inomnot //17 Notional interest rate + mc //19 Real marginal cost + lam //20 Inverse marginal utility of wealth + g //21 Growth shock + s //22 Risk premium shock + mp //23 Monetary policy shock + pi //24 Observed inflation + @#if !(small_model) + x //3 Investment + k //4 Capital + u //7 Utilization cost + ups //8 Utilization choice + wg //9 Real wage growth gap + xg //10 Investment growth + rk //14 Real rental rate + q //18 Tobins q + @#endif + ; + +varexo + epsg // Productivity growth shock + epsi // Notional interest rate shock + epss // Risk premium shock + ; +parameters + // Calibrated Parameters + beta // Discount factor + chi // Labor disutility scale + thetap // Elasticity of subs. between intermediate goods + thetaw // Elasticity of subs. between labor types + nbar // Steady state labor + eta // Inverse frish elasticity of labor supply + delta // Depreciation + alpha // Capital share + gbar // Mean growth rate + pibar // Inflation target + inombar // Steady gross nom interest rate + inomlb // Effective lower bound on gross nominal interest rate + sbar // Average risk premium + // Parameters for DGP and Estimated parameters + varphip // Rotemberg price adjustment cost + varphiw // Rotemberg wage adjustment cost + h // Habit persistence + rhos // Persistence + rhoi // Persistence + sigz // Standard deviation technology + sigs // Standard deviation risk premia + sigi // Standard deviation mon pol + phipi // Inflation responsiveness + phiy // Output responsiveness + nu // Investment adjustment cost + sigups // Utilization + ; + + +// ---------------- Calibration -----------------------------------------// + +beta = 0.9949; // Discount factor +thetap = 6; // Elasticity of subs. between intermediate goods +thetaw = 6; // Elasticity of subs. between labor types +nbar = 1/3; // Steady state labor +eta = 1/3; // Inverse frish elasticity of labor supply +delta = 0.025; // Depreciation +alpha = 0.35; // Capital share +gbar = 1.0034; // Mean growth rate +pibar = 1.0053; // Inflation target +sbar = 1.0058; // Average risk premium +rkbar = 1/beta; +varphiw = 100; // Rotemberg wage adjustment cost +nu = 4; // Investment adjustment cost +sigups = 5; // Utilization +varphip = 100; // Rotemberg price adjustment cost +h = 0.80; // Habit persistence +rhos = 0.80; // Persistence +rhoi = 0.80; // Persistence +sigz = 0.005; // Standard deviation +sigs = 0.005; // Standard deviation +sigi = 0.002; // Standard deviation +phipi = 2.0; // Inflation responsiveness +phiy = 0.5; // Output responsiveness +inomlb = 1 ; // Inom LB + +// ---------------- Model -----------------------------------------------// +model; + +@#if !(small_model) + [name = 'HH FOC utilization (1)'] + rk = steady_state(rk)*exp(sigups*(ups-1)); + + [name = 'Utilization definition (3)'] + u = steady_state(rk)*(exp(sigups*(ups-1))-1)/sigups; + + [name = 'Firm FOC capital (4)'] + rk = mc*alpha*g*yf/(ups*k(-1)); + + [name = 'HH FOC capital (17)'] + q = beta*(lam/lam(+1))*(rk(+1)*ups(+1)-u(+1)+(1-delta)*q(+1))/g(+1); + + [name = 'HH FOC investment (18)'] + 1 = q*(1-nu*(xg-1)^2/2-nu*(xg-1)*xg)+beta*nu*gbar*q(+1)*(lam/lam(+1))*xg(+1)^2*(xg(+1)-1)/g(+1); + + [name = 'Wage Phillips Curve (20)'] + varphiw*(wg-1)*wg = ((1-thetaw)*w+thetaw*wf)*n/yf + beta*varphiw*(lam/lam(+1))*(wg(+1)-1)*wg(+1)*(yf(+1)/yf) ; + + [name = 'Real wage growth gap (6)'] + wg = pigap*g*w/(gbar*w(-1)); + + [name = 'Law of motion for capital (15)'] + k = (1-delta)*(k(-1)/g)+x*(1-nu*(xg-1)^2/2); + + [name = 'Investment growth gap (14)'] + xg = g*x/(gbar*x(-1)); +@#endif + + +@#if small_model + [name = 'Production function (2)'] + yf = n; + + [name = 'Firm FOC labor (5)'] + w = mc*yf/n; + + [name = 'Output definition (7)'] + y = (1-varphip*(pigap-1)^2/2)*yf; + + [name = 'ARC (13)'] + c = y; + + [name = 'Household labpur supply equals flex wage'] + w = wf; +@#else + [name = 'Production function (2)'] + yf = (ups*k(-1)/g)^alpha*n^(1-alpha); + + [name = 'Firm FOC labor (5)'] + w = (1-alpha)*mc*yf/n; + + [name = 'Output definition (7)'] + y = (1-varphip*(pigap-1)^2/2-varphiw*(wg-1)^2/2)*yf - u*k(-1)/g; + + [name = 'ARC (13)'] + x = y-c; +@#endif + +[name = 'Output growth gap (8)'] +yg = g*y/(gbar*y(-1)); + +[name = 'Notional Interest Rate (9)'] +inomnot = inomnot(-1)^rhoi*(inombar*pigap^phipi*yg^phiy)^(1-rhoi)*exp(mp); + +[name = 'Nominal Interest Rate (10)', bind='zlb'] +inom = inomlb; +[name = 'Nominal Interest Rate (10)', relax='zlb'] +inom = inomnot; + +[name = 'Inverse MUC (11)'] +lam = c-h*c(-1)/g; + +[name = 'Flexible real wage definition (12)'] +wf = chi*n^eta*lam; + +[name = 'HH FOC bond (16)'] +1 = beta*(lam/lam(+1))*s*inom/(g(+1)*pibar*pigap(+1)); + +[name = 'Price Phillips Curve (19)'] +varphip*(pigap-1)*pigap = 1-thetap+thetap*mc+beta*varphip*(lam/lam(+1))*(pigap(+1)-1)*pigap(+1)*(yf(+1)/yf); + +[name = 'Stochastic productivity growth (21)'] +g = gbar+ sigz*epsg; + +[name = 'Risk premium shock (22)'] +s = (1-rhos)*sbar + rhos*s(-1)+sigs*epss ; + +[name = 'Notional interest rate shock (23)'] +mp = sigi*epsi; + +[name = 'Observed inflation (24)'] +pi = pigap*pibar; + +end; + +occbin_constraints; +name 'zlb'; bind inom <= inomlb; relax inom > inomlb; +end; + +// ---------------- Steady state -----------------------------------------// +steady_state_model; +mp = 0; +xg = 1; +s = sbar; +n = nbar; +ups =1; +u =0; +q =1; +g = gbar; +yg = 1; +pigap = 1; +wg =1; +pi = pigap*pibar; +// FOC bond +inom = gbar*pibar/(beta*s); +inomnot = inom; +inombar = inom; +// Firm pricing +mc = (thetap-1)/thetap; +// FOC capital +rk = gbar/beta+delta-1; +// Marginal cost definition +@#if small_model +alpha = 0; +thetaw=1; +@#endif +w = (mc*(1-alpha)^(1-alpha)*alpha^alpha/rk^alpha)^(1/(1-alpha)); +@#if small_model +wf = w; +@#else +wf = (thetaw-1)*w/thetaw; +@#endif +// Consolidated FOC firm +k = w*n*gbar*alpha/(rk*(1-alpha)); +// Law of motion for capital +x = (1-(1-delta)/gbar)*k; +// Production function +yf = (k/gbar)^alpha*n^(1-alpha); +// Real GDP +y = yf; +// Aggregate resouce constraint +c = y-x; +// FOC labor +lam = (1-h/gbar)*c; +chi = wf/(n^eta*lam); + +end; + +// ---------------- Checks -----------------------------------------// +steady; +check; + +// ---------------- Simulation -----------------------------------------// +shocks; +var epsi = 1; +var epss = 1; +var epsg = 1; +end; + +steady; +check; + +// ---------------- Estimation -----------------------------------------// + +varobs yg inom pi; +load('dataobsfile','inom') +// check if inom is at lb and remove data + associated shock +verbatim; + inom(inom==1)=NaN; +end; +inx = strmatch('epsi',M_.exo_names); +if any(isnan(inom)) + M_.heteroskedastic_shocks.Qscale_orig.periods=find(isnan(inom)); + M_.heteroskedastic_shocks.Qscale_orig.exo_id=inx; + M_.heteroskedastic_shocks.Qscale_orig.scale=0; +else + options_.heteroskedastic_filter=false; +end + +copyfile dataobsfile.mat dataobsfile2.mat +save dataobsfile2 inom -append + +// -----------------Occbin ----------------------------------------------// +options_.occbin.smoother.debug=1; +occbin_setup(filter_use_relaxation,likelihood_piecewise_kalman_filter); +options_.heteroskedastic_filter=1; +options_.nobs=110; +calib_smoother(datafile=dataobsfile2,first_obs=1,smoothed_state_uncertainty,smoother_redux); +oo_PKF=oo_; + +occbin_setup(likelihood_inversion_filter,smoother_inversion_filter); +calib_smoother(datafile=dataobsfile2,first_obs=1,smoothed_state_uncertainty,smoother_redux); + +figure +subplot(3,1,1) +plot(1:options_.nobs,oo_.SmoothedShocks.epsg,'b-',1:options_.nobs,oo_PKF.SmoothedShocks.epsg) +subplot(3,1,2) +plot(1:options_.nobs,oo_.SmoothedShocks.epsi,'b-',1:options_.nobs,oo_PKF.SmoothedShocks.epsi) +subplot(3,1,3) +plot(1:options_.nobs,oo_.SmoothedShocks.epss,'b-',1:options_.nobs,oo_PKF.SmoothedShocks.epss) +figure +plot(1:options_.nobs,oo_.SmoothedVariables.inom,'b-',1:options_.nobs,oo_PKF.SmoothedVariables.inom) + +burnin=40; +if max(abs(oo_.SmoothedShocks.epsg(burnin+1:end)-oo_PKF.SmoothedShocks.epsg(burnin+1:end)))>0.1 || ... + max(abs(oo_.SmoothedShocks.epsi(burnin+1:end)-oo_PKF.SmoothedShocks.epsi(burnin+1:end)))>0.1 || ... + max(abs(oo_.SmoothedShocks.epss(burnin+1:end)-oo_PKF.SmoothedShocks.epss(burnin+1:end)))>0.1 + error('Smoothed shocks differ too much') +end + +if max(abs(oo_.SmoothedVariables.inom(burnin+1:end)-oo_PKF.SmoothedVariables.inom(burnin+1:end)))>0.001 || ... + max(abs(oo_.SmoothedVariables.yg(burnin+1:end)-oo_PKF.SmoothedVariables.yg(burnin+1:end)))>0.001 || ... + max(abs(oo_.SmoothedVariables.pi(burnin+1:end)-oo_PKF.SmoothedVariables.pi(burnin+1:end)))>0.001 + error('Smoothed observables differ too much') +end + +temp_data=load('dataobsfile2.mat'); +if max(abs(temp_data.inom(options_.first_obs+burnin:options_.first_obs+options_.nobs-1)-oo_PKF.SmoothedVariables.inom(burnin+1:end)))>0.0001 || ... + max(abs(temp_data.yg(options_.first_obs+burnin:options_.first_obs+options_.nobs-1)-oo_PKF.SmoothedVariables.yg(burnin+1:end)))>0.0001 || ... + max(abs(temp_data.pi(options_.first_obs+burnin:options_.first_obs+options_.nobs-1)-oo_PKF.SmoothedVariables.pi(burnin+1:end)))>0.0001 + error('Smoothed observables differ too much') +end +