% Executable 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.
% Data are at a quarterly frequency

%% Upload Dynare
clc ; clear all ;
addpath C:\dynare\6.0\matlab ;

%% Read the data, compute mean and standard deviaton of output growth, and save the output growth data for use in the model code

DataTable = readtable("C:\Work\Input data in quarterly frequency.xls", "Sheet","Data"); % Use the path to your directory

LY = 100 * log(DataTable.Y) ; % Takes the log of output and multiplies by 100 to transform decimal growth into approximate percent change

Dates = datetime(DataTable.observation_date(2:(end),:),"Format","yyyyQQQ") ; % Transforms the observation-date column into a dat-time object

Table_LY_DIF = table( 4 * diff(LY), 'VariableNames', "LY_DIF") ; % Computes q/q growth at annual rate and renames it as in the model code
writetable(Table_LY_DIF, "IN_quarterly_datafile.csv") ; % Saves the data in a csv file, for use in the model code

calc_mean = mean(Table_LY_DIF.LY_DIF,'omitnan') ; % Computes the mean of output growth, ommitting not-available values, for use in the model code as the constant average growth rate g
calc_std = 0.15 * std(Table_LY_DIF.LY_DIF,'omitnan'); % Computes the standard deviation of output growth

%% Run the model

dynare Model_code_quarterly_frequency.mod ; 

%% Accumulate rates of growth to produce the level variables output and trend output

LY_BAR_ACC = cumsum( 0.25 * oo_.SmoothedVariables.LY_BAR_DIF) ;
LY_BAR = LY_BAR_ACC - LY_BAR_ACC(99) + 100 ; % This line sets observation 54 as the base year for trend output, you may set the base year as you wish
LY_HAT = oo_.SmoothedVariables.LY_HAT ;
LY = LY_BAR + oo_.SmoothedVariables.LY_HAT ;
GAMMAA = oo_.SmoothedVariables.GAMMAA ;

%% Accumulate rates of growth to produce the counterfactual variables output and counterfactural trend output

LY_BAR_C_ACC = cumsum( 0.25 * oo_.SmoothedVariables.LY_BAR_DIF_C ) ;
LY_BAR_C = LY_BAR_C_ACC - LY_BAR_C_ACC(99) + LY_BAR(99) ; % This line sets observation 99 as the base year for counterfactual trend output, you may set the base year as you need
LY_BAR_C(1:99-1) = NaN ;

%% Compute output and trend output in deviation form

LY_BAR_D_ACC = cumsum( 0.25 * oo_.SmoothedVariables.LY_BAR_DIF_D ) ;
LY_BAR_D = LY_BAR_D_ACC - LY_BAR_D_ACC(99) ; 
LY_BAR_D(1:99-1) = NaN ;
LY_D = LY_BAR_D + LY_HAT ;

%% Present the filtration results in a figure

figure() ;
subplot(2,2,1)
plot(Dates,[ LY LY_BAR LY_BAR_C]) ;
xlim([Dates(1), Dates(111)]) ;
legend('Output','Trend output','Counterfactual output','Location','Northwest','Box','Off','Fontsize',15) ; 

subplot(2,2,2)
plot(Dates,[ LY_D LY_BAR_D ]) ;
xlim([Dates(98), Dates(111)]) ;
legend('Output in diviation form','Trend output in deviation form','Location','Southeast','Box','Off','Fontsize',15) ; 

subplot(2,2,3)
plot(Dates,[LY_HAT zeros( 111,1 )]) ;
xlim([Dates(1), Dates(111)]) ;
legend('Output gap','Location','Southwest','Box','Off','Fontsize',15) ; 

subplot(2,2,4)
plot(Dates,[GAMMAA g*ones(1,size(GAMMAA,1))' zeros( 111,1 )]) ;
xlim([Dates(1), Dates(111)]) ; % Note that you may need to change the length of the zero vector if you have a different sample
legend('Long-run trend output growth rate','Constant average growth rate','Location','Southwest','Box','Off','Fontsize',15) ; 

%% Save the filtration results in an Excel file

OUT_quarterly = table(LY, LY_BAR, LY_HAT, GAMMAA, LY_BAR_C, LY_D, LY_BAR_D, 'VariableNames', {'Output','Trend output','Output gap','Long-run trend output growth rate','Counterfactual trend output','Output in deviation form','Trend output in deviation form'}) ;
writetable([table(string(Dates), 'VariableNames',"observation_date") OUT_quarterly],"OUT_quarterly.xlsx","WriteVariableNames",true);