#

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