#

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

#

# Written by Dr.
Dimitrios D. Thomakos

#

# e-mail:
thomakos@uop.gr

#

# STSA Version 2.0

#

# Example file #22

#

# Analysis of the
SOI and Recruitment data in Shumway and Stoffer,

# "Time Series
Analysis and its Applications", published by Springer-Verlag.

#

# Initialize

clear

ginit

clc

# Input the data
series

fname1 =
"data\soi.dat"

fname2 =
"data\recruit.dat"

N = 453.

z1 = read(fname1,"double",N,1)

z2 = read(fname2,"double",N,1)

# First,
detrend both series

[b,seb,covb,rss,s2,prd,x]
= trend(z1,1,0)

[b,seb,covb,rss,s2,prd,y]
= trend(z2,1,0)

# Do the analysis of
the recrutiment series

#

# Compute the ACF
and PACF

[ky,ry,se]
= acfplot(y,30,0.05,"Recruitment")

# Compute the
Yule-Walker estimates of an AR(2) model

[phi0,s2] = ar_yw(y,2)

# Use these
estimates as initial values for nonlinear estimation

[phi,e]
= arma_estimate(y,2,0,phi0)

[delta,se_delta,cov_delta,rss,s2,e,test,pval]
= arma_details(y,2,0,phi,20,1)

# Next,
do some forecasting and plot the actuals and forecasts

[yf,msef]
= arma_forecast(y,2,0,phi,s2,24)

# Compute standard
error bands

yfhigh = yf+sqrt(msef)

yflow =
yf-sqrt(msef)

# Merge with last 96
observations and plot

gaddwin("NOTE: the forecasts are for the
detrended series!!!")

gcolor("black")

gxaxis("linear",0,120,10,1)

gtitle("Actual values and 24 forecasts for
recruitment series")

gxgrid("major")

gygrid("major")

gplot({y(end-95::end),yf})

gcolor("red")

gstyle("dashed")

gplot([{y(end-95::end),yflow},{y(end-95::end),yfhigh}])

# Now,
we will estimate a transfer function in the time domain using

# simultaneous
least squares. NOTE: this differs from the published

# approach
in Shumway in Stoffer where they estimate the transfer

# function
parameters one step at a time and it can easily be

# reproduced using
the STSA arma_estimate function

#

# First, we need to
filter both the input and output series

# using an
appropriate ARMA model for the input which is an AR(1)

[phi,s2]
= ar_yw(x,1)

ex = ar_noise(x,phi)

ey = ar_noise(y,phi)

# Next, we compute
and plot the cross-correlation function

[cyx,cxy,se]
= ccf(ex,ey,30,0.05)

plot_ccf(cyx,cxy,se,"REC","SOI")

# Finally
we estimate the parameters using simultaneous nonlinear least squares

# with the model
orders and delay set to the published values...

sv = { 0.1, 0.1, 0.1, 0.1, 0.1 }

[theta,e]
= transfer_estimate(y,x,1,0,5,2,1,sv,0)

print "Estimates of transfer function between
recruitment -y- and soi -x-"

print "Autoregressive estimate for y (omega)
= ",theta(1)

print "Regression estimate for x (delta) =
",theta(2)

print "Autoregressive estimates for noise
phi(1) and phi(2) = ",theta(3::4)'

print "Moving average estimate for noise theta(1) = ",theta(5)

print

print "Note that this differs somehow from
the published solution because of the different method used!"

print "Below we repeat the published
solution..."

print

# Now,
the publishe solution...

#

# Fit
the regression first

yr = y(6::end)

y1 = y(5::end-1)

x1 = x(1::end-5)

[b,seb,covb,rss,s2,R2,R2a,prd,u]
= regress([y1,x1],yr,"First regression of
transfer function modeling");

# Now, apply moving average filter with know coefficient...

eta = ma_noise(u,-0.848)

# Finally, estimate
an AR(2) model for the filtered series eta, again
using the

# above regression function...

yeta = eta(3::end)

xeta = seqlags(eta,2)

regress(xeta,yeta,"Second regression of transfer
function modeling");

# Next, we do the
spectral and co-spectral analysis of the two series

#

# We
compute the spectrum of the SOI series using three different methods...

# First, the
smoothed periodogram estimator using a sine taper

[f1,S1] = spectrum(x,0,10)

plot_spectrum(f1,S1,"linear",1,"SOI - smoothed
periodogram")

# Next, a parametric
autoregressive estimator - the number of lags equals to the

# order
selected by the AICc criterion...

[aic,bic]
= ar_order(x,20,0,1)

paic = aic(1,1)

print "AICc selected lag order = ",paic

[f2,S2] = ar_spectrum(x,paic,0)

plot_spectrum(f2,S2,"linear",1,"SOI - parametric
AR estimator")

# Finally,
a spectrum based on smoothed autocovariances...

[f3,S3] = acvf_spectrum(x,(N-1)/4)

plot_spectrum(f3,S3,"linear",1,"SOI - smoothed
autocovariances")

# Next,
compute the squared-coherence between the two series

[f,Ryx,Syx,Sxy,Sx,Sy]
= squared_coherence(x,y,(N-1)/4,1,"SOI","REC")

# Then,
compute and plot the impulse response coefficients with REC as input!!!

[s,b]
= impulse_response(y,x,40,1,"REC","SOI")