clear
close all
clc

% =========================================================================
% DESCRIPTION 
% This script loads in a FRED-MD dataset, and transform IP, Employment and
% CPI given forecasting horizon for out-of-sample forecast exercise. 
%
% -------------------------------------------------------------------------
% Auxiliary functions
%
%   transform_y() - this subfunction apply log difference transformation 
%   given forecasting horizon. 
%
% -------------------------------------------------------------------------
% INPUT FILES 
%       FRED-MD data in CSV file, e.g. 2015-04.csv or 2022-04.csv. 

% OUTPUT FILES
%       yvar_'year'.mat: MAT file with transformed variables, can be used
%       for out-of-sample forecasts. 
% -------------------------------------------------------------------------

% PARAMETERS TO BE CHANGED

% File name of desired FRED-MD vintage
version = 2024;

if version == 2015
    csv_in = '2015-04.csv';
elseif version == 2024
    csv_in = '2024-12.csv';
end

% specify the forecasting horizons
horizon = [1 6 12];

% =========================================================================
% PART 1: LOAD AND LABEL DATA

% Load data from CSV file
dum=importdata(csv_in,',');

% extract variable names
varnames = dum.textdata(1,2:end);

% Transformation numbers
tcode=dum.data(1,:);

% Raw data
seriesdata=dum.data(2:end,:);

% Dates
dates = dum.textdata(3:end,1);
dates = datetime(dates,'InputFormat','MM/dd/uuuu');

% Month/year of final observation
final_datevec=datevec(dates(end,1));
final_month=final_datevec(2);
final_year=final_datevec(1);

% Month/year of first observation
first_datevec= datevec(dates(1,1));
first_month = first_datevec(2);
first_year = first_datevec(1);

% generate numerical date vector year + (month - 1)/12. e.g. Jan 1963 is
% represented by 1963, and Dec 1963 is 1963+11/12. 
dates_num_yvar = [first_year+(first_month-1)/12:1/12:final_year+(final_month-1)/12]';

% find index of variables to be predicted 
yname = {'INDPRO', 'CPIAUCSL'}; % IP, CPI
[~,Loc_y] = ismember(yname,varnames);

% select out the variables to be predicted
data =  seriesdata(:,Loc_y);

% Transform y variable with function transform_y()

yvar = transform_y(data,horizon);

IP_1 = yvar(:,1,1);
IP_6 = yvar(:,1,2);
IP_12 = yvar(:,1,3);

CPI_1 = yvar(:,2,1);
CPI_6 = yvar(:,2,2);
CPI_12 = yvar(:,2,3);

output_filename = strcat('yvar','_',string(version));
save(output_filename, 'IP_1', 'IP_6', 'IP_12',  'CPI_1', 'CPI_6', 'CPI_12', 'dates', 'dates_num_yvar')

%% Subfunction for Out-of-Sample Inference
% This function apply the appropriate transformation to IP, Employment and
% CPI given forecasting horizon for out-of-sample forecast exercise. 
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% INPUT
    % data:         data to be transformed in a matrix
    % horizon:      all forecasting horizons in a vector
           
% OUTPUT
    % yvar:         transformed variables in a T x N matrix  
    
function yvar = transform_y(data,horizon)

% Dimensions of variables
[a,b] = size(data);                         

L = length(horizon);

% Initialize a temporary matrix to store dependent variables
yvar = NaN*ones(a,b,L);

% Transform data
for j = 1:L
    
    h = horizon(j);
    
    % transform variables by 1200/h*(log(yt+h) - log(yt))
    temp(1+h:a,:) = (1200/h)*(log(data(1+h:a,:)./data(1:a-h,:)));
    
    yvar(:,:,j) = temp;
    
    clear temp
end

end
