#

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

#

# In this example we simulate a nonlinear (exponential) trend with Gaussian MA innovations

# and illustrate how one can use the functions gauss_newton and bhhh_ml with a user-specified

# function for optimization and estimation of a model's parameters.

#

# 1. Plot the realization along with the 'true' exponential trend and a fitted linear

#    quadratic trend for comparison.

# 2. Draw a scatterplot of the series vs. its first lagged value to indicate the strong

#    linearity in adjacent values generated by the exponential trend.

# 3. Use the function gauss_newton to estimate the parameters of the trend function only.

# 4. Analyze the resulting residuals for evidence of remaining serial correlation; as expected,

#    the residuals have serial correlation indicating a MA model.

# 5. Estimate MA parameters of the residuals.

# 6. Then, jointly estimate all the parameters of the model (including the trend parameters,

#    the MA parameters and the innovation variance) using Maximum Likelihood.

#

 

# Clear workspace and command window, initialize graphics, seed the random number generator

clear

clc

ginit

ranseed

 

# Set lenght of realization

N = 200.

 

# Simulate an exponential trend with moving average errors

y = zeros(N,1)           # initialize series

z = zeros(N,1)

x = zeros(N,1)

t = double(seq(int(N))/N)     # create trend variable

s = 0.2                  # std. deviation of Gaussian deviates

e = snormal(N+2,1)*s          # simulate Gaussian deviates

a = 0.5                  # constant parameter for trend

b = 2.5                  # slope parameter for trend

theta = { 0.6, 0.4 }          # moving average coefficients

for i = 3 to N begin

     z(i) = a*exp(b*t(i))

     x(i) = e(i) + e(i-1)*theta(1) + e(i-2)*theta(2)

     y(i) = z(i) + x(i)

end

 

# Crop the series at the beginning

t = t(3::end)

z = z(3::end)

x = x(3::end)

y = y(3::end)

 

# Plot the series, the true exponential trend and a linear quadratic trend for comparison

gaddwin("'True' exp. trend in green and fitted quadr. trend in red")

font_type = "Arial"

font_size = 8

gtitlefont(font_type,font_size)

gxtickfont(font_type,font_size)   

gytickfont(font_type,font_size)

gxtitlefont(font_type,font_size)

gtitle("The data generating process is x(t) = 0.5*exp(2.5*t) + e(t) + 0.6*e(t-1) + 0.4*e(t-2), e(t) ~ N(0,0.2)")

gxaxis("linear",0,N,4)

gxtitle("Time")

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

# Estimate the quadratic trend using the trend function

[b,seb,covb,rss,s2,prd,x,u] = trend(y,2,0)

# Plot the series and the trends

gplot([y,z,prd])

 

# Draw a scatterplot of the y(t) vs. y(t-1)

caption = "Scatterplot"

xtitle = "y(t-1)"

ytitle = "y(t)"

# Note that the scatterplot uses the default kernel and bandwidth values!

scatterplot(y(1::end-1),y(2::end),0,0,caption,xtitle,ytitle)

 

# Estimate the trend coefficients by nonlinear least squares

# using the STSA gauss_newton function

#

# Make y and t global variables

global y

global t

# Define the residuals in a function

function trend_residuals(_gamma) begin

     alpha = _gamma(1)

     _beta = _gamma(2)

     residuals = y - alpha*exp(_beta*t)

return residuals

end

# Set starting values, convergence tolerance and maximum number of iterations

sv    = { 1., 1. }

ht    = 1e-3

tol   = { 1e-4, 1e-3 }

miter = 100.

# Function call; note that this function provides std. errors via the covariance matrix of the estimates!

# and also formatted screen output

#

# Print iteration information and screen output

level = 5.

out   = 1.

[_gamma,S,Cov,iter] = gauss_newton(function trend_residuals,sv,ht,tol,miter,level,out)  

 

# Compute the trend residuals and their variance

u  = trend_residuals(_gamma)

s2 = S/(N-2)

 

# Call the function arma_residual_diagnostics to examine the trend residuals

df = N-2

 m = 30.

 a = 0.05

arma_residual_diagnostics(u,s2,df,m,a)

print

 

# The diagnostics indicate the presence of serial correlation (as expected);

# estimate an ARMA(0,2) model for the residuals

p  = 0.

q  = 2.

sv = { 0.1, 0.1 }

[delta,e] = arma_estimate(u,p,q,sv)

# Get std. errors etc. via the arma_details function and get formatted screen output

out = 1.

[delta,se_delta,cov_delta,rss,s2,e,test,pval] = arma_details(u,p,q,delta,m,out)

print

 

# Now, jointly estimate all the parameters of the process using

# the STSA function bhhh_ml (maximum likelihood)

#

# Construct a function that returns a vector of negative log-likelihoods

function _log_lf(_lambda) begin

     alpha = _lambda(1)

     _beta = _lambda(2)

     u = y - alpha*exp(_beta*t)

     theta = _lambda(3::4)

     sigma = _lambda(5)

     e = ma_noise(u,theta)

     lf = -0.5*log(2*pi) - 0.5*log(sigma^2) - 0.5*(e^2)/(sigma^2)

return -lf

end

 

# Set starting values and call the function; note that we here a larger number of iterations

# as the algorithm can be slow to converge

sv = { 1., 1., 0.1, 0.1, 0.1 }

ht = 1e-3

miter = 1500.

# Ask for formatted output at the end and iteration information every 40 iterations

level = 40.

out   = 1.

[_lambda,S,J,Cov,iter] = bhhh_ml(function _log_lf,sv,ht,tol,miter,level,out)