#

# STATISTICAL TIME SERIES ANALYSIS (STSA) module for O-Matrix.

# Written by Dr. Dimitrios D. Thomakos

# University of Peloponesse

# e-mail: thomakos@uop.gr

# STSA Version 2.0

# Example file #01

# In this example we simulate a stationary, zero mean Gaussian autoregressive model of order 2

# and perform various operations on the resulting realization.

# 1. Plot the realization.

# 2. Estimate and plot the autocorrelation and autocovariance functions.

# 3. Estimate and plot the Fourier spectrum by smoothing the periodogram and using

# and autoregressive approximation.

# 4. Estimate and plot the empirical distribution function of the realization and

# perform two tests for Gaussianity: the Kolmogorov-Smirnov test and the test

# based on the QQ (quantile-quantile) correlation coefficient.

# 5. Estimate the parameters of the 'true' order (=2) and the order of an autoregressive model

# selected by order selection criteria.

# 6. Perform a likelihood ratio test for the statistical significance of the extra estimated

# coefficients.

# 7. Compute the roots of the associated autoregressive polynomial of the model selected

# using the order selection criteria.

# 8. Compute and plot the impulse responses of the model, i.e. the associated MA coefficients

# in the infinite MA representation of the AR model.

# Clear workspace and command window, initialize graphics, seed the random number generator

clear

ginit

clc

# Simulate a realization from a stationary Gaussian ARMA(2,0) process

# Note: you can change the values of the parameters below

mu = 0. # zero mean

phi = { 1.25, -0.75 } # autoregressive coefficients

theta = 0. # moving average coefficients

sigma = 0.2 # innovation std. deviation

N = 500. # length of realization

# Function call, put the realization in the (N * 1) vector x

x = arma_simulate(mu,phi,theta,sigma,N)

# Plot the realization

gaddwin("Sample realization from a Gaussian ARMA(2,0) process")

font_type = "Arial"

font_size = 8

gtitlefont(font_type,font_size)

gxtickfont(font_type,font_size)

gytickfont(font_type,font_size)

gxtitlefont(font_type,font_size)

gtitle("The data generating process is x(t) = 1.25*x(t-1) - 0.75*x(t-2) + e(t), e(t) ~ N(0,0.2)")

gxaxis("linear",0,N,4)

gxtitle("Time")

gcolor({"black","green"})

gplot([x,zeros(N)])

# Plot and save the autocovariance and autocorrelation functions

h = 20 # lag length

alpha = 0.05 # level of significance for std. error bars

varn = "x" # variable name

[kx,rx,se] = acfplot(x,h,alpha,varn)

# Estimate and plot the spectrum

# (a) by smoothing the periodogram using a sine taper

sm = 0. # flag; do not substract the sample mean

K = 20. # number of sine tapers to use

[f1,S1] = spectrum(x,sm,K)

# (b) the spectrum using an autoregressive approximation

p = 2. # "true" autoregressive order; you can change this

sm = 0. # flag; do not substract the sample mean

[f2,S2] = ar_spectrum(x,p,sm)

gaddwin("Fourier spectrum for ARMA(2,0) realization")

gtitle("Green = smoothed periodogram, red = AR approximation")

gxtitle("Frequency")

gxaxis("linear",-0.5,0.5,10)

scale = "log" # select y-axis scale: "linear" or "log"

gyaxis(scale)

gcolor({"green","red"})

gplot(f1,S1)

gplot(f2,S2)

# Examine the Gaussianity of the realization using a QQ plot and the

# Kolmogorov-Smirnov test; test values and pvalues appear in the titles of plots

caption = "QQ Plot for sample realization of a Gaussian ARMA(2,0) process"

plotq = 1. # flag; create the plot

[rQ,sl] = qqplot(x,caption,plotq)

caption = "Empirical and Normal CDF for sample realization of a Gaussian ARMA(2,0) process"

[xvalues,z,w,tests] = Gaussianity_KStest(x,caption,plotq)

# Retrieve the underlying innovation sequence of the process

e = ar_noise(x,phi)

gaddwin("Innovations of sample realization of a Gaussian ARMA(2,0) process")

gplot(e)

# Estimate the coefficients of the realization using the Yule-Walker equations

p = 2.

[phi,s2_H0] = ar_yw(x,p)

print "Yule-Walker coefficient estimates = ", phi'

print "Innovation variance = ", s2_H0

# Estimate the autoregressive model order using order selection criteria

maxp = 10. # Max lag to consider

sm = 0. # flag; do not substract the sample mean

method = 0. # flag; estimation method - use Yule-Walker

[aic,bic] = ar_order(x,maxp,sm,method)

print "AIC order selected = ", aic(1,1)

print "BIC order selected = ", bic(1,1)

print "True order = ", 2

print "Note the PACF plot as well!!!"

# Use the selected AIC order to estimate the coefficients using Yule-Walker

paic = aic(1,1)

[phi,s2_HA] = ar_yw(x,paic)

print "Yule-Walker coefficient estimates using the AIC order = ", phi'

print "Innovation variance = ", s2_HA

# If paic > p then use a Likelihood Ratio test to examine the significance of the

# extra coefficients

if paic > p then begin

# Compute the max of the log-likelihood function under

# the null and under the alternative hypothesis

loglf_H0 = -0.5*N*log(2*pi+1) - 0.5*N*log(s2_H0)

loglf_HA = -0.5*N*log(2*pi+1) - 0.5*N*log(s2_HA)

lr_test = 2*(loglf_HA - loglf_H0)

print "LR test for the significance of the extra coefficients = ", lr_test

print "P-value = ", cdfchi(lr_test,paic-p)

end

# Compute the roots of the autoregressive polynomial of the estimated AIC model

[ar,ma] = arma_roots(paic,0,phi)

print "Estimated abs. values of roots of autoregressive polynomial of AIC model = ", abs(ar)'

print "(Roots should be greater than one for stationarity)"

# Compute and plot the first 30 impulse response coefficients (infinite MA representation) of the estimated AIC model

psi = ar_to_ma(phi,30)

gaddwin("First 30 impulse response coefficients of AIC model")

gtitle("Coefficients in infinite MA representation")

gxtitle("Lag")

gxaxis("linear",0,30,30)

gstyle("diamond")

gcolor("black")

gplot(psi)

gstyle("solid")

gcolor("green")

gplot(psi)