#

# 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 #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")