%% Example of Time Series Analysis

clear all;
close all;
clc;

%% Step 1. Load Time Series.


load('timeseries.mat')
T = length(Y);

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

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

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

% If 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.  The downward sloping of the plot ACF could indicate a unit root process.

figure
subplot(2,2,1)
scatter(Y(2:T),Y(1:T-1))
axis equal
title('lag 1')
subplot(2,2,2)
scatter(Y(3:T),Y(1:T-2))
title('lag 2')
subplot(2,2,3)
scatter(Y(4:T),Y(1:T-3))
axis equal
title('lag 3')
subplot(2,2,4)
scatter(Y(5:T),Y(1:T-4))
title('lag 4')

%% 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] = adftest(Y) 

% THe p-value is 0.001, less than alpha=0.05, so the H0 is rejected with significance level of 5%. The hypotesis of non stationarity of the time series is rejected with level of confidence of 95%.
%% Step 4. Estimate Model.
% Create the model template for estimation.

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

% Here the data are stationary, so we set always D=0, that is the model we use is ARMA.
% 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.
% In statistics, a Q–Q plot (quantile-quantile plot) is a probability plot, a graphical method for comparing two probability distributions by plotting their quantiles against each other. A point (x, y) on the plot corresponds to one of the quantiles of the second distribution (y-coordinate) plotted against the same quantile of the first distribution (x-coordinate). qqplot(x) displays a quantile-quantile plot of the quantiles of the sample data x versus the theoretical quantile values from a normal distribution. If the distribution of x is normal, then the data plot appears linear.
% Function ksdensity(x) returns a probability density estimate, f, for the sample data in the vector x. The estimate is based on a normal kernel function, and is evaluated at equally-spaced points, xi, that cover the range of the data in x.

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 heteroscedasticity and conclude that there are significant ARCH effects in the series.
% Here h=0 and pvalue grater than alpha=0.05, so we do not reject the null hypotesis. The residuals are homoscedastic. Homoscedasticity means that they are roughly the same throughout, which means your residuals do not suddenly get larger.
%
%% Step 9. 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 50 periods.
%

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

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

figure
plot(151:T,y2,'r','LineWidth',2)
hold on
plot(151:T,yF1,'k--','LineWidth',1.5)
xlim([151,T])
title('Prediction Error')
legend('Observed','Forecast','Location','NorthEast')
hold off

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

Mdl(1) = arima(1,0,1);
Mdl(2) = arima(1,0,2);
Mdl(3) = arima(1,0,3);
Mdl(4) = arima(2,0,1);
Mdl(5) = arima(3,0,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(3,0,1) model.