#

#
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 #23

#

# Analysis of the
glacial varve data in Shumway and Stoffer, "Time Series Analysis and its
Applications",

# published by
Springer-Verlag.

#

# Initialize

clear

ginit

clc

# Input the data and
plot the series

fname = "data\varve.dat"

N = 634.

z = read(fname,"double",N,1)

gtitle("Glacial Varve thickness")

gxaxis("linear",0,650,10)

gplot(z)

# Take
log-differences and adjust mean

x = dlog(z)

m = colmean(x)

y = x - m

# Fit an MA(1) model to the log-differenced and mean adjusted series

[delta,e]
= arma_estimate(y,0,1,0.1)

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

print

# As suggested by
Shumway and Stoffer, estimate an ARMA(1,1) model for

# the
log-differenced series as a better model

sv = ones(2,1)*0.1

[delta,e]
= arma_estimate(y,1,1,sv)

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

# Next, we do the
long-memory analysis of the log series

x = log(z)

m = colmean(x)

y = x - m

# Plot the ACF and
PACF of the log series

[kx,rx,se]
= acfplot(x,100,0.05,"Log Glacial Varve")

# Compute the pi
coefficients of the infinite autoregressive representation

[p,z]
= fdiff(x,0.384,0)

gaddwin("Coefficients of infinite autoregressive
representation")

gyaxis("linear",-0.4,0.1,5,2)

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

gcolor("black")

gstyle("square")

gplot(p(2::31));

gstyle("solid")

gplot(p(2::31));

gplot(zeros(30,1));

# Now,
estimate the long memory parameter using nonlinear least squares

# and
the arfima_estimate function

[d,se]
= arfima_estimate(y,0,0,30,0.1,10,1)

print

print "The coefficients in the plot were
generated with the published value of 0.384"

print "which is to be compared with the above
estimate of 0.380..."

print

# Now,
we repeat using frequency domain methods and the Whittle likelihood approach

[d,se_d,loglf]
= fractional_Whittle(y,1,0.2)

print "Estimate
of fractional order by Whittle approximation = ",d,"(",se_d,")"

print "(using FFT)"

print

[d,se_d,loglf]
= fractional_Whittle(y,0,0.2)

print "Estimate
of fractional order by Whittle approximation = ",d,"(",se_d,")"

print "(using DFT)"

print

print
"Published solution = ",0.394,"(0.022)"

print

# Use also the GPH
regression to estimate the fractional order

[d,se_d]
= fractional_GPH(y,0,300)

print "Estimate
of fractional order by GPH regression = ",d,"(",se_d,")"

print