%% ================= Load, Transform (std/mean), Save CSV =================
clearvars; clc; close all;

% --- Paths
thisFileDir  = pwd;
projectRoot  = fileparts(thisFileDir);
dataFolder   = fullfile(projectRoot,'1_Forecasting','data');
outFolder    = fullfile(projectRoot,'1_Forecasting');
if ~isfolder(outFolder), mkdir(outFolder); end

% --- Files to load
fredFiles = { ...
    'INDPRO.xlsx', ...
    'Ricardo_Expectation.xlsx', ...
    'REAINTRATREARAT10Y.xlsx', ...
    'UNRATE.xlsx', ...
    'PCE_NBD20060101.xlsx' ...
};

% --- Load into struct FRED
FRED = struct();
for i = 1:numel(fredFiles)
    filePath = fullfile(dataFolder, fredFiles{i});
    [~, baseName] = fileparts(filePath);

    % Detect sheet (prefer Monthly)
    try, sheets = sheetnames(filePath); 
    catch, [~, sheets] = xlsfinfo(filePath); end
    if isstring(sheets), sheets = cellstr(sheets); end
    preferred = {'Monthly','Monthly Data','Data','Sheet1','Annual'};
    sheetToUse = '';
    for s = 1:numel(preferred)
        if any(strcmpi(sheets, preferred{s})), sheetToUse = preferred{s}; break; end
    end

    % Import
    opts = detectImportOptions(filePath,'Sheet',sheetToUse,...
        'VariableNamingRule','preserve','DatetimeType','datetime');
    opts.VariableNamesRange = 'A1'; opts.DataRange = 'A2';
    T = readtable(filePath, opts);

    % Normalize table
    if ~isdatetime(T{:,1})
        try, datesCol = datetime(T{:,1}, 'ConvertFrom','excel');
        catch, datesCol = datetime(string(T{:,1})); end
    else, datesCol = T{:,1}; end
    valsCol = double(T{:,2});
    FRED.(matlab.lang.makeValidName(baseName)) = table(datesCol, valsCol, ...
        'VariableNames', {'Date','Value'});
end

%% --- Compute std/mean of k-month diffs over 2010s
stdStart = datetime(2010,1,1); stdEnd = datetime(2019,12,1);
K = 1:24;  % horizons
sigmaK = struct(); meanK = struct();

namesLog = {'INDPRO','PCE_NBD20060101'}; % log-transform series
namesRaw = {'Ricardo_Expectation','REAINTRATREARAT10Y','UNRATE'};

for sName = [namesLog, namesRaw]
    T = sortrows(FRED.(sName{1}),'Date');
    maskWin = (T.Date >= stdStart) & (T.Date <= stdEnd);
    valsWin = T.Value(maskWin);
    if ismember(sName,namesLog), valsWin = log(valsWin); end

    sigmaK.(sName{1}) = nan(1,numel(K));
    meanK.(sName{1})  = nan(1,numel(K));
    for k = K
        if numel(valsWin)>k
            diffk = valsWin(1+k:end) - valsWin(1:end-k);
            sigmaK.(sName{1})(k) = std(diffk,'omitnan');
            meanK.(sName{1})(k)  = mean(diffk,'omitnan');
        end
    end
end

%% --- Normalize 2021 relative to Jan 2020
refDate = datetime(2020,1,1);
startDate = datetime(2021,1,1); endDate = datetime(2021,12,31);

% Helper to normalize one series
normalizeSeries = @(vals, refVal, meanK, sigmaK, dates) arrayfun(@(tIdx) ...
    ((vals(tIdx)-refVal-meanK(tIdx)) ./ sigmaK(tIdx)), ...
    1:numel(dates));

% Dates for 2021
dates2021 = FRED.INDPRO.Date(FRED.INDPRO.Date>=startDate & FRED.INDPRO.Date<=endDate);

% Compute k = months since Jan 2020
kIdx = (year(dates2021)-year(refDate))*12 + (month(dates2021)-month(refDate));

% Extract values (log for INDPRO,PCE)
vals.INDPRO = log(FRED.INDPRO.Value(ismember(FRED.INDPRO.Date,dates2021)));
vals.PCE    = log(FRED.PCE_NBD20060101.Value(ismember(FRED.PCE_NBD20060101.Date,dates2021)));
vals.Ricardo_Expectation   = FRED.Ricardo_Expectation.Value(ismember(FRED.Ricardo_Expectation.Date,dates2021));
vals.REAL10 = FRED.REAINTRATREARAT10Y.Value(ismember(FRED.REAINTRATREARAT10Y.Date,dates2021));
vals.UNRATE = FRED.UNRATE.Value(ismember(FRED.UNRATE.Date,dates2021));

ref.INDPRO = log(FRED.INDPRO.Value(FRED.INDPRO.Date==refDate));
ref.PCE    = log(FRED.PCE_NBD20060101.Value(FRED.PCE_NBD20060101.Date==refDate));
ref.Ricardo_Expectation   = FRED.Ricardo_Expectation.Value(FRED.Ricardo_Expectation.Date==refDate);
ref.REAL10 = FRED.REAINTRATREARAT10Y.Value(FRED.REAINTRATREARAT10Y.Date==refDate);
ref.UNRATE = FRED.UNRATE.Value(FRED.UNRATE.Date==refDate);

% Normalize (skip if k outside 1..24)
Y_INDPRO = arrayfun(@(i) (vals.INDPRO(i)-ref.INDPRO-meanK.INDPRO(kIdx(i))) ./ sigmaK.INDPRO(kIdx(i)), 1:numel(kIdx))';
Y_PCE    = arrayfun(@(i) (vals.PCE(i)-ref.PCE-meanK.PCE_NBD20060101(kIdx(i))) ./ sigmaK.PCE_NBD20060101(kIdx(i)), 1:numel(kIdx))';
Y_Ricardo_Expectation   = arrayfun(@(i) (vals.Ricardo_Expectation(i)-ref.Ricardo_Expectation) ./ sigmaK.Ricardo_Expectation(kIdx(i)), 1:numel(kIdx))';
Y_REAL10 = arrayfun(@(i) (vals.REAL10(i)-ref.REAL10) ./ sigmaK.REAINTRATREARAT10Y(kIdx(i)), 1:numel(kIdx))';
Y_UNRATE = arrayfun(@(i) (vals.UNRATE(i)-ref.UNRATE) ./ sigmaK.UNRATE(kIdx(i)), 1:numel(kIdx))';

%% --- Build timetable and export CSV
TT = timetable(dates2021(:), Y_INDPRO(:), Y_Ricardo_Expectation(:), Y_REAL10(:), Y_UNRATE(:), Y_PCE(:));
TT.Properties.VariableNames = {'RealActivity','ExpectedInfl1Y','RealRate10Y','OneMinusUnemp','PCE'};

outFileCSV = fullfile(dataFolder, 'VF_ALL_refJan2020.csv');
writetimetable(TT, outFileCSV);

fprintf('↳ Saved normalized data to %s\n', outFileCSV);
