Introduction to Big Data Analytics

• Week 0: Coding and MATLAB refresher
• Week 1: Monte Carlo Simulation
• Week 2a: Bootstraping
• Week 2b: Cross-Validation (CV): Recursive, Rolling, K-fold
• Week 3: Time Series Decomposition
• Week 4: Outlier Detection and Correction
• Week 5: Penalized Regressions
• Week 6: PCA and PLS

Coding and MATALB Refresher

• Kevin Sheppard's MFE Toolbox: https://www.kevinsheppard.com/code/matlab/mfe-toolbox/
• Interface to download FRED Data: https://www.mathworks.com/matlabcentral/fileexchange/116955-interface-to-easily-access-fred-data

• X-13ARIMA-SEATS Toolbox by Yvan Lengwiler: https://uk.mathworks.com/matlabcentral/fileexchange/49120-x-13-toolbox-for-seasonal-filtering, https://www.census.gov/data/software/x13as.Win_X-13.html
• MIDAS Toolbox by Hang Qian: https://mathworks.com/matlabcentral/fileexchange/45150-midas-matlab-toolbox
• Lasso and elastic-net regularized generalized linear models Toolbox by Hastie, Friedman, Tibshirani and co-authors: https://hastie.su.domains/glmnet_matlab/
• MATLAB Central File Exchange (- repository for 3rd party MATLAB functions): https://mathworks.com/matlabcentral/fileexchange/

Fundamental MATLAB Classes

Week 1: Monte Carlo (MC) Simulation

1. Simulate (i.e. draw) R independent and identically distributed (iid) random samples for X, where X has a know distribution, say . Therefore i = 1, .., R are the number of replications.
2. For every random sample i, evaluate the statistic of interest (e.g. sample mean or variance) .
3. The collection of all provides the approximation to the distribution of the statistic.

clear

%Draw random numbers from a known distribution

T=500;

%1) Standard Normal x~N(0,1)

x = randn(T,1); %Tx1 vector of random scalars drawn from the standard normal distribution

histogram(x, 'Normalization','pdf')

%2) Normal x~N(mu,sigma)

mu=0; sigma=25;

x = mu + sigma*randn(T,1);

%x = normrnd(mu,sigma,T,1) %Same as above

histogram(x, 'Normalization','pdf')

histfit(x) %fits a normal density function together with the histogram

Example 1: Simulate Normal variable - Construct the Mean () Distribution

clear

load('data.mat','y');

plot(y)

xlim([1 length(y)]); title('S&P500 Daily Log-Returns');

T=size(y,1);

R=1000; %# of MC samples

mu=mean(y);

sigma=std(y);

x = mu + sigma*randn(T,R); %TxR matrix of random scalars drawn from the N(mu, sigma) distribution

MCmu = mean(x); %Distribution for the mean estimator

histfit(MCmu,[],'normal')

prctile(MCmu, [2.5 97.5]) %Confidence interval (for mu_hat) at the 95% confidence level

ans = 1×2

-0.0004 0.0831

clearvars -except y

figure();

histfit(y,[],'normal'); title('S&P500 Daily Log-Returns');

%NB: Earlier we plotted the simulated means, here we plot yt itself

Example 2: Simulate AR(1) process - Construct the Distribution

We now want to estimate the vector and then (MC) simulate the process y, multiple times, in order to estimate multiple AR(1) models, one for each of the simulated paths of y, and consequently obtain a distribution of 's.

%------------------------%

% Estimation on actual y %

%------------------------%

%Estimate the parameters of the ARIMA(p,D,q) model via maximum likelihood given the observed time series y

%Notes:

% AR(1) = ARIMA(1,0,0)

% If P ~ ARIMA(1,1,0) then y=diff(P) ~ ARIMA(1,0,0)

mdl = arima('Constant',NaN,'ARLags',1);

%mdl = arima('Constant',NaN,'ARLags',1,'MALags',1) %Also adding MA

[ESTreg, vcov_hat, logLikelihood, info] = estimate(mdl, y);

ARIMA(1,0,0) Model (Gaussian Distribution): Value StandardError TStatistic PValue ________ _____________ __________ ________ Constant 0.046713 0.022536 2.0728 0.038192 AR{1} -0.14964 0.0071289 -20.991 0 Variance 1.1512 0.012025 95.729 0

summarize(ESTreg); %summary of ML regression results

%Extract the estimated coefficients

phi0 = ESTreg.Constant;

phi1 = ESTreg.AR{1};

%ma1 = ESTreg.MA{1}; %If we also had MA component

sigma2 = ESTreg.Variance;

summarize(ESTreg).BIC;

summarize(ESTreg).AIC;

% Now the objective is to SIMULATE y (using Monte Carlo) %

T=size(y,1);

R=1000; %# of MC samples

%-------------------%

% Simulation of y's %

%-------------------%

%Simulate R yt paths, each having T observations

%NB: To be able to simulate, the model-object needs to have the variance specified

mdl = arima('Constant',phi0,'AR',phi1,'Variance',sigma2);

[YR, ER, V] = simulate(mdl,T,'NumPaths',R); % Simulate 1000 y~AR(1)

%Outputs:

%- YR: Simulated paths for y

%- ER: innovations from the ARIMA model

%- V: conditional variances

%-----------------------------%

% Estimation on simulated y's %

%-----------------------------%

%For each of the simulated y paths/columns in YR, estimate the AR(1) model

%and obtain the coefficients

phi = nan(R,2);

for i = 1:R

mdl = arima('Constant',NaN,'ARLags',1) ;

[ESTreg, ~, ~, ~] = estimate(mdl,YR(:,i),'Display','off');

sphi0 = ESTreg.Constant;

sphi1 = ESTreg.AR{1};

phi(i,:) = [sphi0 sphi1];

end

figure()

subplot(2,1,1);

histfit(phi(:,1),[],'normal')

title('Phi0 Distribution')

subplot(2,1,2);

histfit(phi(:,2),[],'normal')

title('Phi1 Distribution')

CI = table(); %table preallocation

CI.Pct = string([2.5; 97.5; 50]);

tmp = prctile(phi, [2.5 97.5 50]); %Confidence interval at the 95% confidence level

CI =[CI array2table(tmp,'VariableNames',["phi0"; "phi1"])];

%Add estimated simulated-mean

CI{end+1,1} = "Mean";

CI{end,2:end} = mean(phi);

%Add estimated phi's [Conditional means: E(phi/yt)] from original sample

CI{end+1,1} = "E(phi/yt)";

CI{end,2:end} = [phi0 phi1];

CI %Return the table

CI = 5×3 table

	Pct	phi0	phi1
1	"2.5"	0.0041	-0.1909
2	"97.5"	0.0874	-0.1120
3	"50"	0.0479	-0.1509
4	"Mean"	0.0472	-0.1505
5	"E(phi/yt)"	0.0467	-0.1496

Potentially the AR(1) assumption for the daily S&P500 log-returns is not a good one, because as we saw from the time-series plot earlier, the actual S&P500 daily returns exhibit volatility clustering. You can repeat the simulation exercise but instead of simulating an AR(1) process (for y), you can now simulate a GARCH(1,1) process and construct the Distribution (of course, now, the 's will be the GARCH coefficients). The GARCH(1,1) model is given by:

with and

.

The μ parameter in the conditional mean equation is called the 'mean-offset'. Modify the code given below to (MC) simulate R paths of as a GARCH(1,1) process. Recall that GARCH(0,1) = ARCH(1).

load('data.mat','y')

figure();

plot(y) %observe the volatility clustering

title('S&P500 Daily Log-Returns');

T = length(y);

xlim([1 T]);

R=1000; %# of MC samples

%------------------------%

% Estimation on actual y %

%------------------------%

mdl = garch('Constant',NaN,'GARCHLags',1,'ARCHLags',1,'Offset',NaN);

[ESTmdl, vcov_hat, logLikelihood, info] = estimate(mdl, y);

GARCH(1,1) Conditional Variance Model with Offset (Gaussian Distribution): Value StandardError TStatistic PValue ________ _____________ __________ __________ Constant 0.046132 0.0044393 10.392 0 GARCH{1} 0.74141 0.016811 44.102 0 ARCH{1} 0.2197 0.016364 13.426 0 Offset 0.07973 0.013702 5.8188 5.9258e-09

phi0 = ESTmdl.Constant;

phi1 = ESTmdl.GARCH{1};

phi2 = ESTmdl.ARCH{1};

mu = ESTmdl.Offset; %parameter 'mu' in the conditional mean equation

summarize(ESTmdl).BIC;

summarize(ESTmdl).AIC;

%-------------------%

% Simulation of y's %

%-------------------%

%Simulate R yt paths, each having T observations

mdl = garch('Constant',phi0,'GARCH',phi1,'ARCH',phi2,'Offset',mu);

[VR, ER] = simulate(mdl,T,'NumPaths',R); % Simulate 1000 y~GARCH(1,1)

%Outputs:

%- VR: Simulated paths for conditional variance (NOT y)

%- ER: Simulated response y. Specifically:

% 1) If (conditional mean) model contains offset (mu), then:

% ER = innovations plus an offset

% 2) If (conditional mean) model does NOT contain offset (mu), then:

% ER = innovations

%NB: You can also simulate by directly passing the estimated model object

%(ESTmdl), instead of specifying it.

[VR, ER] = simulate(ESTmdl,T,'NumPaths',R);

%NB: We simulated R=1000 paths, but let's plot only 2 (out of the 1000).

figure

subplot(2,1,1)

plot(VR(:,1:2))

xlim([0,T])

title('Conditional Variances')

subplot(2,1,2)

plot(ER(:,1:2))

xlim([0,T])

title('Responses (y)')

Week 2: Bootstraping & Cross-Validation (CV)

Bootstraping

1. Draw B random samples with replacement from data X (each of which having the same number of obs as the initial matrix/series X, say T), and collect their bootstrap counterparts {}. Therefore, b = 1, .., B are the number of Bootstrap resamples.
2. For every resample b, evaluate the statistic of interest (e.g. sample mean or variance) .
3. The collection of all provides the bootstrap approximation to the distribution of the statistic.

• The sample autocorrelation function (ACF) and partial autocorrelation function (PACF) are useful qualitative tools to visually assess the presence of autocorrelation at individual lags.
• The Ljung-Box Q-test is a more quantitative way to test for autocorrelation at multiple lags, jointly.

Using the Ljung-Box Q-test you can test for instance the null hypothesis that the first 10 autocorrelations are jointly zero. This is done in MATLAB using function lbqtest() by setting the Number of Lags to . The usual choice for the lags is 20, unless the series has less observations, in which case you can set the lags to be tested equal to T-1. The degrees of freedom (DOF) of the null distribution must be set equal to the number of lags (unless the residuals come from a model with K parameters, in which case, the DOF should be set equal to the number of lags being tested minus K (). The Ljung–Box Q test is appropriate for assessing autocorrelation in any series with a constant mean, therefore, the series needs to be transformed to stationary before passed to the lbqtest(). The Q test statistic (under the null) is defined as follows:

where m denotes the number of autocorrelations calculated (i.e. the number of lags being tested), is the estimated sample autocorrelation for lag j, and T the sample size.

y = fred_api("SP500", 'StartDate',datetime(2005,1,1));

p = rmmissing(y);

p = p{:,1}; %S&P500 daily index

y = diff(log(p))*100; %log-returns

subplot(2,1,1); plot(p); xlim([0,length(p)]); title('S&P500 Index')

subplot(2,1,2); plot(y); xlim([0,length(y)]); title('Daily S&P500 log-return')

%load('data.mat','y')

residuals = y - mean(y); %Compute the centered return series.

T = size(y,1);

lags2test = min([20,T-1]);

[h,pValue,teststat,critValue] = lbqtest(residuals,'DOF',lags2test,'Lags',lags2test,'Alpha',0.05);

%h=1: We reject the no autocorrelation null hypothesis in favor of the alternative.

pValue %print the pValue for the Ljung–Box Q test on daily S&P500 log-returns

pValue = 0

Assume AR(1) process - Construct the Distribution

clear

load('data.mat','y')

%Estimate the AR(1) model using yt

mdl = arima('Constant',NaN,'ARLags',1);

ESTreg = estimate(mdl, y);

summarize(ESTreg); %summary of ML regression results

ARIMA(1,0,0) Model (Gaussian Distribution) Effective Sample Size: 2514 Number of Estimated Parameters: 3 LogLikelihood: -3744.16 AIC: 7494.33 BIC: 7511.82 Value StandardError TStatistic PValue ________ _____________ __________ ________ Constant 0.046713 0.022536 2.0728 0.038192 AR{1} -0.14964 0.0071289 -20.991 0 Variance 1.1512 0.012025 95.729 0

%Extract the estimated coefficients

phi0 = ESTreg.Constant;

phi1 = ESTreg.AR{1};

%------------------------------------------%

% Create the distributions for phi0 & phi1 %

%------------------------------------------%

T=size(y,1);

B=1000; %# of Bootstrap samples

w = ceil(T^1/3);

b8t = block_bootstrap((1:T)',B,w); %MFE Toolbox function

phi = nan(B,2);

for i = 1:B %loop over the B bootstrap samples

disp(['Bootstrap Iter i=', num2str(i), ' out of B=', num2str(B)]);

smpl = b8t(:,i);

mdl = arima('Constant',NaN,'ARLags',1); %yt ~ AR(1)

ESTreg = estimate(mdl,y(smpl),'Display','off');

b_phi0 = ESTreg.Constant;

b_phi1 = ESTreg.AR{1};

phi(i,:) = [b_phi0 b_phi1];

end

figure

subplot(2,1,1);

histfit(phi(:,1),[],'normal')

title('Phi0 Distribution')

subplot(2,1,2);

histfit(phi(:,2),[],'normal')

title('Phi1 Distribution')

CI = table(); %table preallocation

CI.Pct = string([2.5; 97.5; 50]);

tmp = prctile(phi, [2.5 97.5 50]); %Confidence interval at the 95% confidence level

CI =[CI array2table(tmp,'VariableNames',["phi0"; "phi1"])];

%Add estimated phi mean

CI{end+1,1} = "Mean";

CI{end,2:end} = mean(phi);

%Add estimated phi's [Conditional means: E(phi/yt)] from original sample

CI{end+1,1} = "E(phi/yt)";

CI{end,2:end} = [phi0 phi1];

CI %Return the table

CI = 5×3 table

	Pct	phi0	phi1
1	"2.5"	0.0327	-0.2506
2	"97.5"	0.0626	0.0028
3	"50"	0.0456	-0.1498
4	"Mean"	0.0463	-0.1315
5	"E(phi/yt)"	0.0467	-0.1496

Cross-Validation (CV): Recursive, Rolling, K-fold

1. Recursive CV

2. K-fold CV

1. Split the data in X (and y) into K=10 groups, each of which has approximately the same number of observations.
2. Use the 1st group (out of the K groups) of data as the test set (Ztest = ytest,Xtest). Use the remaining (K-1) 9 groups of data as training set (Ztrain = ytrain,Xtrain).
3. Compute the loss [e.g. (R)MSE, MAE] over Ztest using the model estimated over Ztrain.
4. Repeat steps 2 & 3 but with the rest groups. E.g. use the 2nd group of data as the Ztest. Use the remaining as Ztrain.
5. Proceed in a similar manner until each group is used as Ztest set exactly once.
6. Return the 10 computed values (MSE's MAE's etc) as the vector values.

The way to go about optimizing hyperparameters is usually by creating a grid of potential values and and then CV'ing over each of those values.

• For each of the K test/train splits of the sample:
• (1) loop over all lambda values , and estimate the LASSO (obtained using the training part of the sample).
• (2) Predict the Y for the TEST set (using the betas estimated over the TRAINING set).
• Repeat the procedure for each of the K sets.
• You will then obtain a matrix of CV MSEs (ntests x L) (where ntests=K).
• Averaging column-wise (i.e. taking the mean of K values) you will obtain the so-called CROSS-VALIDATION ERROR. It will be 1-by-L vector. Essentially, the CV Error will be the Average of the K Mean MSE's for each λ.
• We select the (optimal) lambda by picking the one with the minimum CV-Error.

Forecasting Performance Evaluation

To compare candidate models, we need to create a measure reflecting how well each of the those competing models perform in an out-of-sample (OOS) forecasting setup. We refer to those as 'error measures'. There are multiple such measures. Denote the pseudo-out-of-sample (POOS) forecast error (FE) coming from model m, with .

%------------------------------------------%

%Simulation of popular stochastic processes%

%------------------------------------------%

T=1000;

R=100; %# of MC samples

Ys = nan(T,1);

%%% WN = ARMA(0,0) %%%

%y(t) = a + e(t) -> ARMA(0,0)

mdl = arima('Constant',0.5,'Variance',1);

%mdl = arima(0,0,0); %Equivalent to above

[YR, ER, V] = simulate(mdl,T,'NumPaths',R);

Ys(:,1) = YR(:,1);

%Outputs:

%- YR: Simulated paths for y

%- ER: innovations from the ARIMA model

%- V: conditional variances

%%% AR(1) with b>0 %%%

%y(t) = a + b*y(t-1) + e(t) -> AR(1)

mdl = arima('Constant',0.5,'Variance',1,'AR',0.8);

[YR, ER, V] = simulate(mdl,T,'NumPaths',R);

Ys(:,2) = YR(:,2);

%%% AR(1) with b<0 %%%

%y(t) = a + b*y(t-1) + e(t) -> AR(1)

mdl = arima('Constant',0.5,'Variance',1,'AR',-0.8);

[YR, ER, V] = simulate(mdl,T,'NumPaths',R);

Ys(:,3) = YR(:,3);

%%% MA(1) %%%

%y(t) = a + e(t) + h*e(t-1) -> MA(1)

mdl = arima('Constant',0.5,'Variance',1,'MA',0.2);

[YR, ER, V] = simulate(mdl,T,'NumPaths',R);

Ys(:,4) = YR(:,4);

%%% ARMA(1,1) %%%

%y(t) = a + b*y(t-1) + e(t) + h*e(t-1) -> ARMA(1,1)

mdl = arima('Constant',0.5,'Variance',1,'AR',0.8,'MA',0.2);

[YR, ER, V] = simulate(mdl,T,'NumPaths',R);

Ys(:,5) = YR(:,5);

%%% RW %%%

Mdl_RW = arima('D',1,'Constant',0,'Variance',1);

[YR, ER, V] = simulate(Mdl_RW,T, 'Y0',0,'NumPaths',R);

Ys(:,6) = YR(:,6);

%%% GARCH(1,1) with offset (offset = constant 'a' in the conditional-mean equation) %%%

%y(t) = a + e(t) -> ARMA(0,0)

%e(t) = sigma(t)Z(t) with Zt~N(0,1)

%sigma(t)^2 = k + g*sigma(t-1)^2 + g*e(t-1)^2 -> GARCH(1,1)

CMmdl = arima('Constant',2.5); %ARMA(0,0) conditional mean model

CVmdl = garch('Constant',0.05,'GARCH',0.8,'ARCH',0.15); %GARCH(1,1) conditional variance model

%Recall: For stationarity we want: (GARCH + ARCH < 1)

CMmdl.Variance = CVmdl;

[YR, ER, V] = simulate(CMmdl,T,'NumPaths',R); %Simulate R y~GARCH(1,1) with offset

Ys(:,7) = YR(:,7);

figure

M = 7;

subplot(M,1,1); plot(Ys(:,1)); xlim([0,T]); title('WN')

subplot(M,1,2); plot(Ys(:,2)); xlim([0,T]); title('AR(1) b=0.8')

subplot(M,1,3); plot(Ys(:,3)); xlim([0,T]); title('AR(1) b=-0.8')

subplot(M,1,4); plot(Ys(:,4)); xlim([0,T]); title('MA(1)')

subplot(M,1,5); plot(Ys(:,5)); xlim([0,T]); title('ARMA(1,1)')

subplot(M,1,6); plot(Ys(:,6)); xlim([0,T]); title('RW')

subplot(M,1,7); plot(Ys(:,7)); xlim([0,T]); title('GARCH(1,1)')

Week 3: Time Series Decomposition (Anomaly Detection I)

In the main lecture we decomposed the time series into all its main components. Specifically, a time series, , can be have the following 3 components:

• Trend component, Tt
• Seasonal component, St with some known periodicity s
• Irregular (stationary) stochastic component, It

.

addpath(genpath(fullfile(pwd, 'MyFunctions'))); %Add folders to MATLAB search-path

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Download SA & NSA series from FRED %

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%1. Industrial Production (Monthly)

TT = fred_api(["INDPRO","IPB50001N"], 'StartDate', datetime(2006,1,1));

TT.Properties.VariableNames = ["SA","NSA"];

figure();

p = stackedplot(TT);

p.Title = 'Industrial Production';

plot(TT.Time, TT{:,1}, TT.Time, TT{:,2})

title('Industrial Production'); legend(["SA","NSA"])

%2. Vehicle Miles Traveled (Units: Millions of Miles; Monthly)

TT = fred_api(["VMTD11","TRFVOLUSM227NFWA"], 'StartDate', datetime(2006,1,1));

TT.Properties.VariableNames = ["SA","NSA"];

p = stackedplot(TT);

p.Title = 'Vehicle Miles Traveled';

plot(TT.Time, TT{:,1}, TT.Time, TT{:,2})

title('Vehicle Miles Traveled'); legend(["SA","NSA"])

%-> Peak in July's; Trough in February's

tr = timerange(datetime(2018,1,1),datetime(2019,12,31));

TT = TT(tr,:);

plot(TT.Time, TT{:,1}, TT.Time, TT{:,2})

title('Vehicle Miles Traveled - 2-year Window'); ; legend(["SA","NSA"])

%3. Quarterly Series: Real GDP

%Real Gross Domestic Product: Billions of Chained 2012 Dollars, SA Annual Rate (GDPC1) -> 2020Q2 $17282.1

%Real Gross Domestic Product: Billions of Chained 2012 Dollars, NSA (ND000334Q) -> 2020Q2: $4329.3

[TT,Meta] = fred_api(["GDPC1","ND000334Q"], 'StartDate', datetime(2006,1,1));

TT.Properties.VariableNames = ["SA","NSA"];

p = stackedplot(TT);

p.Title = 'US Real GDP';

%plot(TT.Time, TT{:,1}, TT.Time, TT{:,2})

• STL (Seasonal and Trend Decomposition using Loess) filtering,
• X13-ARIMA seasonal adjustment program of the U.S. Census Bureau.

Under follows a stationary process.

→ The true process is a RW (i.e. is a unit root process = contains stochastic trend)

As such, the process under :

Different software use different versions of the the ADF regression. In some cases you will find the LHS of the ADF regression to be instead of . However the test result remains the same.

• When regressing on , test .
• When regressing on , test .

• 14400 for monthly data
• 1600 for quarterly data
• 100 for yearly data

TT = fred_api("GDPC1", 'StartDate',datetime(2010,1,1), 'EndDate',datetime(2019,12,30));

[trend,cyclical] = hpfilter(TT{:,1},1600);

TC = array2timetable([trend cyclical],"RowTimes",TT.Time,"VariableNames",["Trend","Cyclical"]);

p = stackedplot([TT TC],{["GDPC1" "Trend"] "Cyclical"});

p.Title = 'HP-filter on US Real GDP';

p.DisplayLabels = {'GDP','Cycle'};

p.AxesProperties(1).LegendLabels = ["GDP","Trend"];

p.AxesProperties(1).LegendLocation = 'SouthEast';

TT = fred_api("GDPC1");

[trend,cyclical] = hpfilter(TT{:,1},1600);

TC = array2timetable([trend cyclical],"RowTimes",TT.Time,"VariableNames",["Trend","Cyclical"]);

plot(TT.Time, TT{:,1}, TT.Time, TC{:,1});

title("HP-filter on US Real GDP"); legend(["GDP","Trend"]); axis tight

% Changing the Smoothing Parameter %

TT = fred_api("GDPC1", 'StartDate',datetime(2008,1,1));

smoothing = [1 50 100 1600 14400 Inf]; %grid for smoothing parameter lambda

tiledlayout(2,3)

for j = 1:numel(smoothing)

nexttile

[trend,cyclical] = hpfilter(TT{:,1},smoothing(j));

TC = array2timetable([trend cyclical],"RowTimes",TT.Time,"VariableNames",["Trend","Cyclical"]);

title("\lambda = " + string(smoothing(j)));

axis tight

plot(TT.Time, TT{:,1}, TT.Time, TC{:,1});

end

clear

Week 4: Outlier Detection and Correction (Anomaly Detection II)

%Create noisy data with 3 outliers:

N=150;

z = 15*rand(N,1);

x = sin(z) + 0.5*(rand(N,1)-0.5);

x(ceil(N*rand(2,1))) = 3; %Two outilers (-the high value of the outlier is 3)

x(ceil(N*rand(1,1))) = 4; %1 more outlier (-taking the value of 4)

figure();

scatter(z,x);

histogram(x, 'Normalization','pdf');

%-----------------------%

%%% Outlier Detection %%%

%-----------------------%

%A. Outliers defined as any value exceeding 2 St.Dev's from the median

%This method will ensure removal of the extreme 5% ONLY if normality can be

%assumed for the series.

%Recall that for a normal population:

%- 68% of the population is within 1 std of the mean.

%- 95% of the population is within 2 std of the mean.

%- 99.7% of the population is within 3 std of the mean.

w = 2;

toobig = median(x) + w*std(x);

toosml = median(x) - w*std(x);

outlier = (x >= toobig) | (x <= toosml); %Indicator/logical vector for the detected outliers

%B. Outliers defined as observations that deviate from the sample

% median by more than 10 interquartile (IQR) ranges

iqr = prctile(x,75) - prctile(x,25); %Interquartile range (IQR). prctile() omits NaNs by default

xmm = abs( x - median(x,'omitnan') );

thr = 10; %threshold for outliers detection (i.e. multiples of IQR)

outlier = (xmm > thr*iqr);

%------------------------%

%%% Outlier Adjustment %%%

%------------------------%

%1. In case you only have a single series x

x_adj = x;

x_adj(outlier) = median(x,'omitnan');

[outlier x_adj x];

%2. If the predictor-matrix X is of size T-by-N.

for i=1:N

x = X(:,i);

if sum(outlier)==0 %no outliers detected

continue

else %outliers detected

if adjm==0 %set outliers to missing

x(outlier) = NaN;

elseif adjm==1 %Median

x(outlier) = median(x,'omitnan');

end

end

X(:,i) = x

end

Week 5: Variable Selection & Variable Reduction: Penalized Regression

where , nestsing LASSO () and RR ().

Therefore, EN is a hybrid between the two abovementioned approaches, the Ridge Regression, and the LASSO regularization. The difference between LASSO and RR is that LASSO coefficients can be equal to 0 (with more 0's for larger λ values), but RR does not set coefficients exactly to 0 (regardless of the λ value), however, coefficients become very-very small. Note also that for small λ values the LASSO coefficients are close to the least-squares (OLS) estimates . Empirical studies suggest that the EN can outperform LASSO on data with highly correlated predictors.

The FRED-MD Monthly Macroeconomic Dataset

• tcode=1 → x~I(0); No transformation: x(t)
• tcode=2 → x~I(1); First difference: x(t)-x(t-1)
• tcode=3 → x~I(2); Second difference: (x(t)-x(t-1))-(x(t-1)-x(t-2))
• tcode=4 → lnx~I(0); Natural log: ln(x)
• tcode=5 → lnx~I(1); First difference of natural log: ln(x)-ln(x-1)
• tcode=6 → lnx~I(2); Second difference of natural log: (ln(x)-ln(x-1))-(ln(x-1)-ln(x-2))
• tcode=7 → x~I(2); First difference of percent change: (x(t)/x(t-1)-1)-(x(t-1)/x(t-2)-1)

Forecasting Inflation

• idx1 = minimum CV-MSE,
• idx2 = minimum CV-MSE plus 1 standard error,

then, by construction, (idx2) will be larger than (idx1) . This implies that idx2 corresponds to:

1. a higher CV-MSE,
2. sparsest model (i.e. less X's picked).

%%% PLOT 1: CV-MSE VS (decending) lambda

%Plot the cross-validated (CV) MSE for the different lambda values.

figure()

lassoPlot(B,S,'PlotType','CV','PredictorNames',Xnames); %-> plots CV MSE for each lambda

legend('show')

%%% PLOT 2: estimated betas VS (decending) lambda

%Plot the non-zero coefficients in the regression for various values of the Lambda

figure()

lassoPlot(B,S,'PlotType','Lambda','XScale','log','PredictorNames',Xnames);

legend('show', 'Location','bestoutside')

By default lasso() will center and scale (i.e. standardize) the predictors for , and will center (i.e. de-mean) the response . However, even though the LASSO betas are estimated using the standardized X's, most functions (including MATLAB's lasso() and glmnet()) will return the estimated LASSO betas on the original data scale. Keep in mind that if you come across any regularized regression function that returns the scaled betas instead, you should then standardize your X's or un-scale the estimated betas, before calculating your predictions.

Can we improve our forecasting model(s)?

ML Interpretability

Week 6: Variable Selection & Variable Reduction: PCA and PLS

Principal Component Analysis (PCA)

The k PCA factors of matrix X are simply the first k eigenvectors of the matrix , with Z denoting the normalized data matrix X. Alternatively, PCA factors and factor loadings can be computed by extracting the first k eigenvectors of the covariance matrix of the normalized dataset X. Denoting the i-th normalized indicator with for , and with Z, the matrix that contains all the normalized indicators, then, the covariance matrix is given by .

Applying PCA on a matrix X with dimensions , will result in:

• T factors in total, if , or
• N factors in total, if .

This implies that the F matrix will be of dimensions or . However, our objective is to reduce the dimension of the problem be selecting the first k. One method for selecting the optimal number of principal components (), is by setting the desired level of the timeseries' variance explained. Alternatively, information criteria (ICs) can be used, such as the generalizations of Mallow's proposed by Bai & Ng (2002), which has over 4k citations.

However, note that those methods for selecting the optimal number of factors () do not take into account the response variable y. Once has been established, (say for instance ) we can then start involving the response variable y and conduct a cross-validation (CV; recursive or rolling) exercise or use information criteria in a regression-based framework (such as F-AR(p)) in order to choose which of the factors (and possibly which of the p lags of those factors) better forecast our variable of interest. Depending on how the CV exercise is set up, this can result in multiple factors being excluded. In other words, the pervious step only tells us the number of factors () we should consider, meaning that we do not necessarily have to include all factors into our final model. For example, the recursive-CV exercise or BIC might conclude that only the 1-st and 5-th factors out of the factors need to be included in our model when forecasting y. The factor-augmented regression (F-AR) framework takes the form:

where and .

A Real-Time Monthly Dataset

%Download a set of Monthly & Daily series via the FRED API

addpath(genpath(fullfile(pwd, 'MyFunctions'))); %Add folders to MATLAB search-path

%The 3rd-party MFE toolbox contains a function named pca().

%We need to remove it from the search-path in order to use MATLAB's built-in pca() function.

rmpath(fullfile('.\MyFunctions\mfe-toolbox-main\crosssection'))

%25 US Monthly Macro Series

%-> Downloaded in real-time using the FRED API

fred_codes = ["CONSUMER";"RPI";"DPCERA3M086SBEA";"INDPRO";"IPMANSICS";"TCU";"CUMFNS";"PAYEMS";"SRVPRD";"UNRATE";"HOUST";"CMRMTSPL";"PCEPI";"PPICMM";"CPIAUCSL";"CPILFESL";"MCOILWTICO";"FEDFUNDS";"TB3MS";"GS10";"TWEXAFEGSMTH";"EXUSUK";"EXUSEU";"UMCSENT";"USEPUINDXM"];

[TT, Meta] = fred_api(fred_codes, 'StartDate',datetime(1995,1,1));

%Take a look on the selected indicators

Meta.Title;

%Let's compare Headline and Core CPI inflation

%NB: Core CPI excludes food and energy prices

infl = ["CPIAUCSL"; "CPILFESL"];

y = 12*100*diff(log(TT{:,infl})); %Annualized MoM% change

TTinfl = array2timetable(y, 'RowTimes',TT.Time(2:end), 'VariableNames',infl);

figure()

plot(TTinfl.Time, TTinfl{:,:});

title('US CPI Inflation (Annualized MoM%)');

legend(["Headline","Core"]); axis tight

p = stackedplot(TTinfl);

p.Title = 'US CPI Inflation (Annualized MoM%)';

p.DisplayLabels = ["Headline","Core"];

clear fred_codes infl y TTinfl p

%Keep the date we downloaded the data

downdatefredAPI = char(datetime('now', 'Format','yyyyMMdd'));

%Saving a row (unprocessed) version of the data is always a good idea when using APIs!

save(fullfile('./Data/fred_raw.mat'));

Xt = TT;

dates = Xt.Time; %Extract dates

%Select the target variable (yt) to forecast

Yname = "CPIAUCSL"; %Headline CPI

%Transform the LHS

y_untr_api = Xt(:,Yname); %Keep a copy of the untransformed LHS

y = nan(size(Xt,1),1);

y(2:end) = 12*100*diff(log(Xt{:,Yname})); %Annualized MoM% change

TTy = array2timetable(y(2:end), 'RowTimes',Xt.Time(2:end), 'VariableNames',Yname);

TTy(end-4:end,1) %Display the last 5 obs

ans = 5×1 timetable

	Time	CPIAUCSL
1	01/05/2022	11.6291
2	01/06/2022	15.7630
3	01/07/2022	-0.2316
4	01/08/2022	1.4175
5	01/09/2022	4.6227

%Xt(:, Yname) = []; %NB: For PCA we do NOT have to remove yt from X matrix!

clear Yidx TT TTy

X = Xt{:,:}; %Convert timetable to array

[T,N] = size(X);

%--------------------------%

% Treating the Ragged-edge %

%--------------------------%

%Unlike in the FRED-MD, there might be NaNs for 'CPIAUCSL' at the end of the sample

%depending on when we are downloading the data (because some of the other indicators

%are released with less delay than CPI). The missing values at the end of

%the sample are called 'ragged-edge' or 'jagged-edge'.

LHS_Tidx = find(~isnan(Xt{:,Yname}), 1,'Last');

LHS_T = Xt.Time(LHS_Tidx);

raggededge = T - LHS_Tidx; %How many NaNs there are at the end of my LHS

clear T LHS_Tidx

%Adjust/re-align the matricies to account for the ragged edge.

%Here I decided to retain the ragged-edge (i.e. the last few rows of X),

%because there is potentially important timely information in those last obs.

%Alternatively, I could have removed those last few rows containing NaNs.

%NB: We will now be forecasting yT+1 using XT+1 (instead of XT), or even using XT+2.

%NB: The code works even if there is no ragged-edge in 'CPIAUCSL', as it

% will set raggededge to 0, and therefore, no re-alignment will take place.

y = y(1:end-raggededge);

X = X(1+raggededge:end,:);

dates = dates(1:end-raggededge);

[T,N] = size(X);

%-----------------%

% Transformations %

%-----------------%

%Which variables to log before transforming?

%Recall: log(-3)=complex; log(0)=Inf; log(1)=0; xt/0=Inf

% => if x<=0 do NOT log()

%SW(2002) log all non-negative series that were not already in rates or percentage units.

%Rule 1: Series that are not negative or 0

small = 1e-6; % Value close to zero

canlog = (min(X) > small); %series NOT close to 0

%Rule 2: Series that are not already in rates or percentage units

%Simply go through the table to see which variables are already in %

Meta(:, ["SeriesID","Title","Units"])

ans = 25×3 table

	SeriesID	Title	Units
1	"CONSUMER"	"Consumer Loans, All Commercial Banks"	"Billions of U.S. Dollars"
2	"RPI"	"Real Personal Income"	"Billions of Chained 2012 Dollars"
3	"DPCERA3M086SBEA"	"Real personal consumption expenditures (chain-type quantity index)"	"Index 2012=100"
4	"INDPRO"	"Industrial Production: Total Index"	"Index 2017=100"
5	"IPMANSICS"	"Industrial Production: Manufacturing (SIC)"	"Index 2017=100"
6	"TCU"	"Capacity Utilization: Total Index"	"Percent of Capacity"
7	"CUMFNS"	"Capacity Utilization: Manufacturing (SIC)"	"Percent of Capacity"
8	"PAYEMS"	"All Employees, Total Nonfarm"	"Thousands of Persons"
9	"SRVPRD"	"All Employees, Service-Providing"	"Thousands of Persons"
10	"UNRATE"	"Unemployment Rate"	"Percent"
11	"HOUST"	"New Privately-Owned Housing Units Started: Total Units"	"Thousands of Units"
12	"CMRMTSPL"	"Real Manufacturing and Trade Industries Sales"	"Millions of Chained 2012 Dollars"
13	"PCEPI"	"Personal Consumption Expenditures: Chain-type Price Index"	"Index 2012=100"
14	"PPICMM"	"Producer Price Index by Commodity: Metals and Metal Products: Primary Nonferrous Metals"	"Index 1982=100"
⋮

pos = [6 7 10 18 19 20]; %Manually add the indicies

idx = false(N,1); %Convert into logical vector

idx(pos) = true;

%Alternatively, use string fuctnion contains() to quantify this decision/selection

idx = contains(Meta.Units, "percent",'IgnoreCase',true);

notpct = ~idx'; %tilda gives the negation operation -> changes 1 to 0, and 0 to 1

%Combine rules 1 & 2

takelog = notpct & canlog; %same as: notpct.*canlog

[notpct; canlog; takelog]; %Check the decisions and which X's will be logged.

X(:,takelog) = log(X(:,takelog)); %Take the natural logarithm for the selected series

clear small canlog pos idx notpct takelog

%Transform to Stationary based on repeating Unit Root (UR) tests

%(since we do not have any transformation codes given to us)

%NB: We dont know how many times we will need to difference each

%of the series until stationarity is achieved, therefore, we need

%to use a WHILE loop instead of FOR loop.

Xst = nan(size(X));

countdiffs = zeros(1,N);

for i=1:N

isStationary = 0;

x = X(:,i);

while isStationary == 0

[h,~,~,~,ADFreg] = adftest(x, 'Model','ARD', 'Lags',0:5, 'Alpha',0.1);

allBIC = [ADFreg.BIC]; %The BIC's

[~, idx] = min(allBIC);

%Keep the results of the ADF model with the BIC-selected lag order

isStationary = h(idx); %h=1 implies stationarity at the 10% level of significance

if isStationary == 0

x = diff(x); %First differences

countdiffs(i) = countdiffs(i) + 1; %Count how many times we have differenced the series

end

end

Xst(1+countdiffs(i):end, i) = x;

end

countdiffs; %Take a look on how many differences it took to stationar-ize each indicator

clear h ADFreg allBIC idx i countdiffs x isStationary

X = Xst; %Make X the transformed matrix

clear Xst

%------------------------------%

% Outlier Detection/Adjustment %

%------------------------------%

iqr = prctile(X,75) - prctile(X,25); %Interquartile range (IQR). prctile() omits NaNs by default

xmm = abs( X - median(X,'omitnan') );

thr = 10; %threshold for outliers detection (i.e. multiples of IQR)

outlier = (xmm > thr*iqr);

clear iqr xmm thr

sum(outlier); %column-wise sum to check # of problematic obs per indicators

%Step 1: Replace detected outliers with missing

for i=1:N

out_i = outlier(:,i);

if sum(out_i)==0 %no outliers detected

continue

else %outliers detected

X(out_i,i) = NaN; %Replace with missing

end

end

clear out_i outlier

%Before we fix the outliers we could remove the rows with too many NaNs

row_out = sum(isnan(X),2); %row-wise summation of units

table(dates, row_out); %NB: We will have a spike in the number of outliers/NaNs around the COVID-lockdowns period

idx = 2; %Remove the upper-edge up to the idx-th observation (inclusive)

X = X(idx+1:end,:);

y = y(idx+1:end);

dates = dates(idx+1:end);

[T,N] = size(X);

clear idx row_out

%Step 2: Correct Outliers and the NaNs

for i=1:N

tofill = isnan(X(:,i));

X(tofill,i) = median(X(:,i),'omitnan');

end

clear i tofill T N ans

%Now the dataset will be 'balanced' (meaning with no missing values).

%NB: Here we fixed the NaNs by using the median. Given the ragged-edge

% problem, and because we are estimating a FAR forecasting model,

% a better fix would be to impute the missing using EM ala SW(2002).

%Save the dataset for future use

save(fullfile('./Data/fred.mat')) %saves the entire workspace

clear

Principal Component Analysis (PCA)

clear

load(fullfile('./Data/fred.mat'));

%SYNTAX:

%[coeff,score,latent,tsquar,explained,mu] = pca(X);

%If X is TxN (with T>N), then:

%- coeff (NxN): Principal component coefficients, also known as loadings

%- score (TxN): Principal component scores are the representations of X in the

% principal component space (i.e. the factors)

%- latent (Nx1): Latent are the principal component variances (i.e. the eigenvalues of the covariance matrix of X)

%- tsquar (Tx1): Hotelling's T-squared statistic for each observation in X

%- explained (Nx1): Percentage of the total variance explained by each principal component

%- mu (1xN): The estimated mean of each variable in X.

%NB: By default pca() centers X by subtracting the means of each column before

% computing the principal components. To apply PCA on the non-centered data, use:

% [coeff,score,latent,tsquared,explained,mu] = pca(Xf, 'Centered', false);

%%% Standardized X (i.e. instead of centered) %%%

%Because MATLAB only de-means (i.e. centers) the X matrix, I can standardize X before applying pca().

Z = (X - mean(X))./std(X);

[coeff,F,~,~,explained] = pca(Z, 'Centered',false);

%Factors can be thought of as weighted-averages of the individual series

%Let's compare the 1st PCA factor with the cross-sectional average using equal weights

T = size(X,1);

plot(1:T, F(:,1), 1:T, mean(Z,2));

title("PCA factor Vs Averaged series"); axis tight

legend(["1st PCA factor", "Time-series Average (of Z)"], 'location','SW');

%Re-construct the (centered) X matrix (i.e. Zhat):

Xcentered = F*coeff'; %-> Zhat

%Let's take a look of yt vis-a-vis the first few factors

TTy = array2timetable(y, 'RowTimes',dates, 'VariableNames',Yname);

Fnames = "F"+ string(1:size(F,2));

TTf = array2timetable(F, 'RowTimes',dates, 'VariableNames',Fnames);

TT = [TTy TTf];

p = stackedplot(TT,{["F1" "F2"] Yname});

p.Title = 'Factors & CPI Inflation (Annualized MoM%)';

p.DisplayLabels = {'F','Infl.'};

p.AxesProperties(1).LegendLabels = ["F1","F2"];

p.AxesProperties(1).LegendLocation = 'SouthWest';

clear TTy Fnames TTf TT p

%-----------------------------------------------%

%%% Variability Explained by PCA and k Choice %%%

%-----------------------------------------------%

explained; %<- Displays the % of total variance explained by the principal components

%Table with Cummulative Variance Explained by each of the first 10 factors

kk = 10;

array2table([(1:kk)' cumsum(explained(1:kk))], 'VariableNames', ["F", "Cummul. Var. Expl."])

ans = 10×2 table

	F	Cummul. Var. Expl.
1	1	17.5947
2	2	31.6041
3	3	41.9913
4	4	50.8909
5	5	57.5251
6	6	63.6013
7	7	68.5373
8	8	72.7965
9	9	76.3515
10	10	79.8102

clear kk

%Make a plot of the incremental percent variability explained by each principal component.

figure()

pareto(explained)

xlabel('Principal Component')

ylabel('Variance Explained (%)')

%Find the number of factors (k) required to explain at least 95% variability.

sum_explained = 0;

k = 0;

while sum_explained < 95

k = k + 1;

sum_explained = sum_explained + explained(k);

end

k %You could use this figure as your k_star (-number of factors to consider in your FAR)

k = 17

%Alternativly, find the k that minimizes the PCp2 information criterion by Bai-Ng (2002)

kmax = 8; %Select the max number of factor to consider

Zhat = Xcentered;

e = Z - Zhat; %Residuals from predicting Z using the factors

sigma = mean(sum(e.*e/T)); %Sum of squared residuals divided by NT

kz = 1:kmax;

PCp2 = log(sigma) + (N+T/N*T)*log(min([N T]))*kz; %Criterion PC_p2

[minPCp2, kstr] = min(PCp2,[],2);

clear Xcentered explained sum_explained kmax Zhat e sigma kz PCp2 minPCp2

Forecasting with PCA Factors