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