clear
clc

% Purpose:
%
% This is a demo of GARCH-MIDAS model
%
% Reference:
%
% Engle,R.F., Ghysels, E. and Sohn, B. (2013), Stock Market Volatility 
% and Macroeconomic Fundamentals. The Review of Economics and Statistics,
% 95(3), 776-797.
%
% Written by Hang Qian
% Contact: matlabist@gmail.com

% NASDAQ Composite Index, daily percentage change 1971 - 2015
% Data Source: FRED Economic Data
% https://research.stlouisfed.org/fred2/series/NASDAQCOM
y = xlsread('SP500.xls','B2343:B11669');

% Estimate the GARCH-MIDAS model, and extract the volatilities
period = 22;
numLags = 24;
[estParams,EstParamCov,Variance,LongRunVar] = GarchMidas(y,'Period',period,'NumLags',numLags); %#ok<*ASGLU>

% Plot the conditional volatility and its long-run component
figure(1)
nobs = size(y,1);
seq = (period*numLags+1:nobs)';
year = linspace(1971.1,2015.8,nobs);
plot(year(seq),sqrt(252*Variance(seq)),'g--','LineWidth',1);
hold on
plot(year(seq),sqrt(252*LongRunVar(seq)),'b-','LineWidth',2);
legend('Total Volatility','Secular Volatility','Location','SouthEast')
xlim([1980,2016])
ylim([0,0.5])
hold off

% Estimate the rolling window version of the GARCH-MIDAS model
[estParams,EstParamCov,Variance,LongRunVar]...
    = GarchMidas(y,'Period',period,'NumLags',numLags,'RollWindow',1);

% Plot the conditional volatility and its long-run component
figure(2)
nobs = size(y,1);
seq = (period*numLags+1:nobs)';
year = linspace(1971.1,2015.8,nobs);
plot(year(seq),sqrt(252*Variance(seq)),'g--','LineWidth',1);
hold on
plot(year(seq),sqrt(252*LongRunVar(seq)),'b-','LineWidth',2);
legend('Total Volatility','Secular Volatility','Location','SouthEast')
xlim([1980,2016])
ylim([0,0.5])
hold off

% Unemployment rate, 1980-2015
% Data Source: FRED database 
% https://research.stlouisfed.org/fred2/series/INDPRO
xMonth = xlsread('UNEMP_USA.xls','b149:b576')./100;

% Repeat the monthly value throughout the days in that month
[~,yDate] = xlsread('SP500.xls','A2343:A11669');
[~,yDateMonth] = datevec(yDate);
xDay = NaN(nobs,1);
count = 1;
for t = 1:nobs
    if t > 1 && yDateMonth(t) ~= yDateMonth(t-1)    
        count = count + 1;
        if count > length(xMonth)
            break
        end
    end
    xDay(t) = xMonth(count);
end

% Estimate the rolling window version of the GARCH-MIDAS model
[estParams,EstParamCov,Variance,LongRunVar] = GarchMidas(y,'Period',period,'NumLags',32,'X',xDay,'ThetaM',1);

% Plot the conditional volatility and its long-run component
figure(3)
nobs = size(y,1);
seq = (period*numLags+1:nobs)';
year = linspace(1971.1,2015.8,nobs);
plot(year(seq),sqrt(252*Variance(seq)),'g--','LineWidth',1);
hold on
plot(year(seq),sqrt(252*LongRunVar(seq)),'b-','LineWidth',2);
legend('Total Volatility','Secular Volatility','Location','SouthEast')
xlim([1980,2016])
ylim([0,0.5])
hold off