#
#
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 #09
#
# In
this example we illustrate the differences between the two available kernels
for nonparametric
# estimation
and forecasting. Then, we illustrate the differences in the fit between a
linear and a
# nonparametric
regression when the true model is nonlinear but unknown.
#
# Clear workspace
and command window, initialize graphics, seed the
random number generator
clear
clc
ginit
ranseed
# Plot the optimal
and Gaussian kernels
x = -3.5:0.1:3.5 # grid of x values
x0 = 1. #
evaluation point; you can change this
h =
0.1 # bandwidth; you can change
this
z =
(x-x0)/h # standardized variate
y2 = gaussian_kernel(z) # the
Gaussian kernel
y1 = optimal_kernel(z) # the
optimal kernel
# Now
plot them
gaddwin("The optimal (green) and Gaussian (red)
kernel functions")
format double "f4.2"
title = ["Evaluation point is x0 =
",ntoa(x0)," and bandwidth is h = ",ntoa(h)]
gtitle(title)
gcolor({"green","red"})
gxaxis("linear",0.,1.5)
gplot(x,[y1,y2])
# Now,
repeat for a different bandwidth value
h =
1. # bandwidth;
you can change this
z =
(x-x0)/h # standardized variate
y2 = gaussian_kernel(z) # the
Gaussian kernel
y1 = optimal_kernel(z) # the
optimal kernel
gaddwin("The kernels with a 10-times larger
bandwidth")
format double "f4.2"
title = ["Evaluation point is x0 =
",ntoa(x0)," and bandwidth is h = ",ntoa(h)]
gtitle(title)
gcolor({"green","red"})
gxaxis("linear",-2.5,4.)
gplot(x,[y1,y2])
# Simulate
a nonlinear regression model and estimate it using
# nonparametric
regression
N = 200. # length
of realization
s = 0.2 # standard deviation of error
variates
e = snormal(N,1)*s #
the regression error variates
x = snormal(N,1) #
the explanatory variates
a = 1. # parameter;
you can change this
b = 5. # parameter;
you can change this
y = a*tanh(-b*x) + e #
the dependent variable
gaddwin("Simulated y = a*tanh(-b*x) + e")
gplot(y)
h = 0.1 # set the bandwidth; you can
change this
# Draw the
scatterplot to show that there is nonlinearity
scatterplot(x,y,0,h,"q","x","y")
# Estimate the
nonparametric regression and a linear one for comparison
z = [ones(N,1),x]
[m,u,stat]
= np_cmean(y,z,0,h) #
nonparametric
[b,seb,covb,rss,s2,prd,u]
= lsq(z,y) #
linear
# Now
plot the observations and the two fits
gaddwin("Observations, linear fit and
nonparametric fit")
gtitle("Circle = observations, red = linear
fit, green = nonparametric fit")
gcolor({"black","red","green"})
gstyle("circle")
gplot(y)
gstyle("solid")
gplot([prd,m])