#

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

#

# In this example we illustrate the estimation and forecasting of ARFIMA-type models.

# This class of models extends traditional, short-memory ARMA models, to incorporate the effects

# of long-memory.

#

# In the example, we illustrate how to use the arfima_estimate and arfima_forecast STSA functions

# and we compare the forecasting performance of ARFIMA forecasts vs. forecasts generated by a linear

# AR model.

#

# Clear workspace and command window, initialize graphics

clear

clc

ginit

ranseed

# Input simulated data

fname = "data\fract.prn"

N = 400.

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

# We can estimate the fractional order and the autoregressive coefficient of the y series

# using an ARFIMA(1,d,0) model...Note that the fractional order of the y series is the

# same as the fractional order of the x series.

#

sv = { 0.1, 0.1 } # Set starting values

# Function call, get the estimates

[theta,S,cov_theta,e,iter] = arfima_estimate(y,1,0,20,sv,5,1)

print

print "True parameter values were = ",[0.5,0.4]

print

# Now, forecast the last few values of the series using the ARFIMA model and an AR model

# and compare their performance...

#

nf = 10. # Number of points to forecast; this can be changed...

yin = y(1::end-nf)

yot = y(end-nf+1::end)

# Estimate and forecast using the AR model

maxp = 10. # max lag to consider; set it high enough as the ACF decays slowly

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

method = 1. # flag; use LS for estimation

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

paic = aic(1,1) # get the selected lag order

sv = ones(paic,1)*0.1

[phi,e] = arma_estimate(yin,paic,0,sv)

s2 = (e'*e)/(N-nf-paic) # get the residual variance

[yf_ar,msef_ar] = arma_forecast(yin,paic,0,phi,s2,nf)

#

# Now estimate and forecast using ARFIMA model

#

sv = { 0.1, 0.1 }

[theta,S,cov_theta,e,iter] = arfima_estimate(yin,1,0,20,sv,5,1)

yf_arf = arfima_forecast(yin,1,0,theta,20,nf)

#

# Compare forecasting performance

#

print

fe_ar = yot - yf_ar

fe_arf= yot - yf_arf

[mae_ar,mse_ar] = foreval(fe_ar)

[mae_arf,mse_arf] = foreval(fe_arf)

print "MAE of AR model = ", mae_ar

print "MAE of ARFIMA model = ", mae_arf

print

print "MSE of AR model = ", mse_ar

print "MSE of ARFIMA model = ", mse_arf

# Now, plot the actual values and the forecasts from the two models

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

gstyle("star")

gtitle("Actuals (stars) and forecasts (red = AR, green = ARFIMA)")

gplot(yot)

gstyle("solid")

gplot(yf_ar)

gplot(yf_arf)