% 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 at annual rate
% DIF4         Four-quarter difference
% Suffix C     Counterfactual
% Suffix D     In deviation from counterfactural
% g            Constant avererage growth rate 
% E_LY_HAT     Demand shock
% 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_BAR_DIF4, LY_BAR_DIF4_D, LY_DIF_D ; 
varexo      E_LY_HAT, E_LY_BAR, E_GAMMAA ;
parameters  alphaa, thetaa, g ;
alphaa = 0.8    ;
thetaa = 0.85   ; % 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 + 4 * E_LY_BAR ;
  GAMMAA      = thetaa * GAMMAA(-1) + ( 1 - thetaa ) * g + E_GAMMAA ;
  LY_DIF      = 4 * ( LY_HAT - LY_HAT(-1) ) + LY_BAR_DIF ;
  LY_BAR_DIF4 = 0.25 * ( LY_BAR_DIF + LY_BAR_DIF(-1) + LY_BAR_DIF(-2) + LY_BAR_DIF(-3) ) ;

  LY_BAR_DIF_C = GAMMAA ;
  
  LY_BAR_DIF_D  = 4 * E_LY_BAR ;
  LY_BAR_DIF4_D = 0.25 * ( LY_BAR_DIF_D + LY_BAR_DIF_D(-1) + LY_BAR_DIF_D(-2) + LY_BAR_DIF_D(-3) ) ;
  LY_DIF_D      = 4 * ( 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_DIF4 = g ;
    LY_BAR_DIF_C = g ;
    LY_BAR_DIF_D = 0 ;
    LY_BAR_DIF4_D = 0 ;
    LY_DIF_D = 0 ;
end ;

shocks;
  var  E_LY_HAT ; stderr dev ;
  var  E_LY_BAR ; stderr 0.42 * dev / 0.8 ;
  var  E_GAMMAA ; stderr 0.25 * dev;
  corr E_LY_HAT, E_LY_BAR = 0.8 ;
end;

varobs LY_DIF ;

%estimated_params;
%   stderr E_LY_BAR, 0.42 * dev / 0.8, inv_gamma_pdf, 0.42 * dev / 0.8, 1.5 ;
%end;
 
%estimated_params_init(use_calibration);
%end;

heteroskedastic_shocks; % Make sure to input your preferred heteroskedastic periods and factors
   var     E_LY_HAT;
   periods 51:58,99:103;
   scales  1.3, 13.7;

   var     E_LY_BAR;
   periods 51:58,99:103;
   scales  1.3, 13.7;
end;

calib_smoother(datafile  = IN_quarterly_datafile,
                heteroskedastic_filter) ;

%estimation(datafile  = IN_quarterly_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_BAR_DIF4, LY_BAR_DIF4_D, LY_DIF_D ;