%% Example of Time Series Analysis

clear all;
close all;
clc;

%% Step 1. Generate Synthetic Time Series.

rng('default') % For reproducibility

% https://it.mathworks.com/help/econ/arima.html?searchHighlight=arima&s_tid=srchtitle_arima_1

constant = 0.2; % Model constant
beta = []; % Regression component coefficients Example: [0.5 NaN 3]
variance = 5; %Model innovations variance Example: 1, Example: garch(1,0)
ar = {0.3}; % nonseasonal AR polynomial coefficients Example: 'AR',{0.5 –0.1},'ARLags',[1 4]
arLags = [1]; % Lags associated with nonseasonal AR polynomial coefficients
ma= {-0.2 -0.5}; % nonseasonal MA polynomial coefficients
maLags = [1 3]; % Lags associated with nonseasonal MA polynomial coefficients
sar= {}; % seasonal AR polynomial coefficients
sarLags = []; % Lags associated with seasonal AR polynomial coefficients
sma= {}; % easonal MA polynomial coefficients
smaLags = []; % Lags associated with seasonal MA polynomial coefficients
d = 1; % degree of nonseasonal integration
seasonality = 0; % degree of seasonal differencing polynomial
m = 0.0; % Trend stationary coefficient

SDMdl = arima('Constant', constant, 'Beta', beta, 'Variance', variance, ...
    'AR', ar,'ARLags', arLags, ...
    'MA',ma,'MALags',maLags, ...
    'SAR', sar, 'SARLags', sarLags, ...
    'SMA',sma,'SMALags',smaLags, ...
    'D', d, 'Seasonality', seasonality);

T = 200;
Y = simulate(SDMdl,T) + m * (1:T)';  

figure;
plot(Y);
title('DATASET')
xlabel('t')
ylabel('Y(t)')

%
%% Step 2. Plot the sample ACF and PACF.
% Plot the sample autocorrelation function (ACF) and partial autocorrelation function (PACF).

figure
subplot(2,1,1)
autocorr(Y)
subplot(2,1,2)
parcorr(Y)

%
% The downward sloping of the plot ACF indicates a unit root process. The lengths of the line segments on the ACF plot gradually decay, and continue this pattern for increasing lags. This behavior indicates a nonstationary series.
%% Step 3. Assess Stationarity Statistically.
% A unit root (also called a unit root process or a difference stationary process) is a stochastic trend in a time series, sometimes called a “random walk with drift”; If a time series has a unit root, it shows a systematic pattern that is unpredictable.Unit root tests are a tool for assessing the presence of a stochastic trend in an observed series.
% Use |adftest| to perform the Dickey-Fuller test on the data. Dickey-Fuller tests assess the null hypothesis of a unit root process in a univariate time series y. Values of h equal to 0 indicate a failure to reject the unit-root process (H0 non stationary). Values of h equal to 1 indicate rejection of the unit-root process in favor of the alternative model (H1 stationary).

display('Dickey-Fuller test for stationarity')
[h,pValue,stat,cValue] = adftest(Y) 

%% Step 4. Estimate Model.
% Create the model template for estimation.

P =1;
D =1;
Q =3;
Mdl = arima(P,D,Q);

% Fit the model to Y.

EstMdl = estimate(Mdl,Y);

%% Step 5. Infer the residuals.

res = infer(EstMdl,Y);
stres = res/sqrt(EstMdl.Variance);

%% Perform Residual Diagnostic Checks
%% Step 6. Check the residuals for autocorrelation.
% Confirm that the residuals are uncorrelated. Look at the sample autocorrelation function (ACF) and partial autocorrelation function (PACF) plots for the standardized residuals.



figure
subplot(2,2,1)
plot(stres)
title('Standardized Residuals')
subplot(2,2,2)
histogram(stres,10)
title('Standardized Residuals')
subplot(2,2,3)
autocorr(stres)
subplot(2,2,4)
parcorr(stres)

display('Ljung-Box Q-test')
[h,p] = lbqtest(stres)

%
% The Ljung-Box Q-test confirms the sample ACF and PACF results.
% The value h = 1 means reject the null hypothesis  of no autocorrelation at the \alpha = 0.05 significance level. The value h = 0 means not reject the null hypothesis  of no autocorrelation at the \alpha = 0.05 significance level (this is white noise).
%
%% Step 7. Check the residuals for normality.
% One assumption of the fitted model is that the residuals follow a Gaussian distribution. Check them for normality.

figure
subplot(1,2,1)
qqplot(stres)

x = -4:.05:4;
[f,xi] = ksdensity(stres);
subplot(1,2,2)
plot(xi,f,'k','LineWidth',2);
hold on
plot(x,normpdf(x),'r--','LineWidth',2)
legend('Standardized Residuals','Standard Normal')
hold off

% The quantile-quantile plot (QQ-plot) and kernel density estimate show no obvious violations of the normality assumption.
%% Step 8. Check the residuals for heteroscedasticity.
% Compute the deviations of the series.

deviations = stres - mean(stres);

%
% Test the series for ARCH effects using the residuals.

display('archtest')
[h,pValue,stat,cValue] = archtest(deviations)

% The result h = 1 indicates that you should reject null hypothesis of no conditional heteroscedasticity and conclude that there are significant ARCH effects in the series.
%
%
%% Check predictive performance.
% Use a holdout sample to compute the predictive MSE of the model. Use the first 100 observations to estimate the model, and then forecast the next 44 periods.

y1 = Y(1:150);
y2 = Y(151:end);

yF1 = forecast(EstMdl,50,'Y0',y1);
pmse = mean((y2-yF1).^2)

figure
plot(y2,'r','LineWidth',2)
hold on
plot(yF1,'k--','LineWidth',1.5)
xlim([0,50])
title('Prediction Error')
legend('Observed','Forecast','Location','NorthWest')
hold off

% The predictive ability of the model is quite good. You can optionally compare the PMSE for this model with the PMSE for a competing model to help with model selection.
%% Compute All Information Criteria
%
% For each competing model, create an arima model template for estimation.

Mdl(1) = arima(1,1,3);
Mdl(2) = arima(2,1,3);
Mdl(3) = arima(5,1,2);
Mdl(4) = arima(3,1,1);

% Fit each model to the simulated data y, compute the loglikelihood, and suppress the estimation display.

numMdl = numel(Mdl);
logL = zeros(numMdl,1);      % Preallocate
numParam = zeros(numMdl,1);

for j = 1:numMdl
    [EstMdl,~,logL(j)] = estimate(Mdl(j),Y,'Display','off');
    results = summarize(EstMdl);
    numParam(j) = results.NumEstimatedParameters;
end

% For each model, compute all available information criteria.

[~,~,ic] = aicbic(logL,numParam,T)

% ic is a 1-D structure array with a field for each information criterion.  Each field contains a vector of measurements; element j corresponds to the model yielding loglikelihood logL(j).
% For each criterion, determine the model that yields the minimum value.

[~,minIdx] = structfun(@min,ic);
[Mdl(minIdx).Description]'

% The minimum of each criterion corresponds to the ARIMA(1,1,3) model
%%
%%