% Model code for 
% Gómez-Pineda, Javier and Roa-Rozo, Julián, A trend-cycle decomposition with hysteresis, Borradores de Economía, No. 1230, March 2023.
%
% LY_HAT       Output gap 
% LY_BAR       Trend output
% GAMMAA       Long-run trend output growth rate
% DIF          First difference 
% Suffix C     Counterfactual
% Suffix D     In deviation from counterfactural
% g            Constant avererage growth rate 
%
% E_LY_HAT     Demand shocks
% E_LY_BAR     Supply-level shock
% E_GAMMAA     Supply-growth shock         
%
% Note that the codes do not use rightly-spelled greek letters to avoid use of possible Matlab function names

var         LY_HAT, GAMMAA, LY_DIF, LY_BAR_DIF, LY_BAR_DIF_C, LY_BAR_DIF_D, LY_DIF_D ; 
varexo      E_LY_HAT, E_LY_BAR, E_GAMMAA ;
parameters  alphaa, thetaa, g, dev ;
alphaa = 0.5 ;
thetaa = 0.9999 ; % Calibrate theta close to 1 in long samples
g      = calc_mean  ;
dev = calc_std ;

model(linear) ;
  LY_HAT = alphaa * LY_HAT(-1) + E_LY_HAT ;
  LY_BAR_DIF = GAMMAA + E_LY_BAR ;
  GAMMAA = thetaa * GAMMAA(-1) + ( 1 - thetaa ) * g + E_GAMMAA ;
  LY_DIF = ( LY_HAT - LY_HAT(-1) ) + LY_BAR_DIF ;

  LY_BAR_DIF_C = GAMMAA ;  
  LY_BAR_DIF_D = E_LY_BAR ;
  LY_DIF_D = ( LY_HAT - LY_HAT(-1) ) + LY_BAR_DIF_D ;
end;

steady_state_model ;
LY_HAT = 0 ;
LY_BAR_DIF = g ;
GAMMAA = g ;
LY_DIF = g ;
LY_BAR_DIF_C = g ;
LY_BAR_DIF_D = 0 ;
LY_DIF_D = 0 ;
end ;

shocks;
  var  E_LY_HAT ; stderr dev ;
  var  E_LY_BAR ; stderr 0.67 * dev / 0.8 ;
  var  E_GAMMAA ; stderr 0.25 * dev ;
  corr E_LY_HAT, E_LY_BAR = 0.7 ;
end;

varobs LY_DIF ;

%estimated_params;
%  stderr E_LY_BAR, 0.67 * dev / 0.8, inv_gamma_pdf, 0.67 * dev / 0.8, 1.5 ;
%end;

%estimated_params_init(use_calibration);
%end;

heteroskedastic_shocks;

   var     E_LY_HAT;
   periods 1:33,57:58,68:69;
   scales  2.2,3.2,5.2;

   var     E_LY_BAR;
   periods 1:33,57:58,68:69;
   scales  2.2,3.2,5.2;

end;

calib_smoother(datafile=IN_annual_datafile,heteroskedastic_filter) ;

%estimation(datafile  = IN_annual_datafile,
%           mode_compute = 5,
%           heteroskedastic_filter,mh_replic=20000) LY_HAT, GAMMAA, LY_DIF, LY_BAR_DIF, LY_BAR_DIF_C, LY_BAR_DIF_D, LY_DIF_D ; 
