GARCHestimation.m

%% GARCHestimation.m

% Supporting script for my blog post about GARCH models.


%   Petr Javorik (2016) maple@mmquant.net


%   http://mmquant.net/introduction-to-volatility-models-with-matlab-garch


%% input data
C = BFX_day1_OHLCV(:,4);
date = BFX_day1_date;

%% Returns. Note that we don't know return for C(1) so we drop first element
r = double((log(C(2:end)./C(1:end-1)))*100); % scaled returns in [%] for numerical stability
e = r - mean(r); % innovations after simple linear regression of returns
C = C(2:end);
date = date(2:end);

%% Plot C and r
% C
figure1 = figure;
subplot1 = subplot(2,1,1,'Parent',figure1);
hold(subplot1,'on');
plot(date,C,'Parent',subplot1);
ylabel('Closing price');
box(subplot1,'on');
set(subplot1,'FontSize',16,'XMinorGrid','on','XTickLabelRotation',45,'YMinorGrid','on');
% r
subplot2 = subplot(2,1,2,'Parent',figure1);
hold(subplot2,'on');
plot(date,r,'Parent',subplot2);
ylabel('returns [%]');
box(subplot2,'on');
set(subplot2,'FontSize',16,'XMinorGrid','on','XTickLabelRotation',45,'YMinorGrid','on');

%% Autocorrelation of returns innovations - ACF, PACF, Ljung-Box test
% ACF
figure2 = figure;
subplot3 = subplot(2,1,1,'Parent',figure2);
hold(subplot3,'on');
autocorr(e); % input to ACF are innovations after simple linear regression of returns
% PACF
subplot4 = subplot(2,1,2,'Parent',figure2);
hold(subplot4,'on');
parcorr(e); % input to ACF are innovations after simple linear regression of returns
% Ljung-Box test
[hLB,pLB] = lbqtest(e,'Lags',3);

%% Conditional heteroskedasticity of returns - ACF, PACF, Engle's ARCH test
% ACF
figure3 = figure;
subplot5 = subplot(2,1,1,'Parent',figure3);
hold(subplot5,'on');
autocorr(e.^2);
% PACF
subplot6 = subplot(2,1,2,'Parent',figure3);
hold(subplot6,'on');
parcorr(e.^2);
% ARCH test
[hARCH,pARCH] = archtest(e,'lags',2);

%% AR-GARCH model, ARIMA object
MdlG = arima('ARLags',2,'Variance',garch(1,1)); % normal innovations
MdlT = arima('ARLags',2,'Variance',garch(1,1)); % t-distributed innovations
MdlT.Distribution = 't';

%% Parameters estimation
% normal innovations
EstMdlG = estimate(MdlG,r);
% t-distributed innovations
EstMdlT = estimate(MdlT,r);

%% Volatility inference and log-likelihood objective function value from estimated AR-GARCH model
[~,vG,logLG] = infer(EstMdlG,r);
[~,vT,logLT] = infer(EstMdlT,r);

%% Fitted models comparison using AIC, BIC
% AIC,BIC
% inputs: values of loglikelihood objective functions for particular model, number of parameters
% and length of time series
[aic,bic] = aicbic([logLG,logLT],[5,6],length(r))

%% AR-GJR-GARCH, ARIMA object
MdlGJR_T = arima('ARLags',2,'Variance',gjr(1,1));
MdlGJR_T.Distribution = 't';

%% Parameters estimation
% t-distributed innovations
EstMdlGJR_T = estimate(MdlGJR_T,r);

%% Volatility inference from estimated AR-GJR-GARCH model
[~,v_GJR_T,logL_GJR_T] = infer(EstMdlGJR_T,r);

%% Fitted models comparison using BIC, AIC
[aic2,bic2] = aicbic([logLT,logL_GJR_T],[6,7],length(r));

%% plot results
% Closing prices
figure4 = figure;
subplot7 = subplot(2,1,1,'Parent',figure4);
hold(subplot7,'on');
plot(date,C);
ylabel('Closing price');
set(subplot7,'FontSize',16,'XMinorGrid','on','XTickLabelRotation',45,'YMinorGrid','on','ZMinorGrid',...
    'on');
% volatility AR-GARCH, innovations t-distributed
subplot8 = subplot(2,1,2,'Parent',figure4);
hold(subplot8,'on');
plot(date,vT);
% volatility AR-GARCH, innovations normally distributed
plot(date,vG);
ylabel('volatility');
legend({'$\varepsilon_t$ $t$-distributed','$\varepsilon_t$ normally distributed'},'Interpreter','latex');
set(subplot8,'FontSize',16,'XMinorGrid','on','XTickLabelRotation',45,'YMinorGrid','on','ZMinorGrid',...
    'on');