
AR1SIM.OMS
Script File:
# Description:
# Simulates an AR(1) process and estimates its spectrum.
#
# An AR(1) process satisfies the recursive formula
# x(i) = a * x(i  1) + w(i)
# where w(n) is white noise. We compare the function's true spectrum
# to a spectral estimator. References below are to
# "Spectral Analysis and Time Series" by M.B. Priestley.
#
clear
# the AR coefficient
a = .8
# the covariance at lags 0 and 1 (see Equation 3.5.16)
r = {1., a} / (1  a^2)
# number of realizations of the process
M = 10
# number of points in each realization
N = 100
# compute the realizations
x = arcov(N, M, r)
# fractional bandwidth for the DPSS taper
w = 2. / N
# compute the dpss taper
taper = dpss(N, 1, w)
# time between data points
dt = 1
# compute the spectral estimate
estimate = fspec(x, dt, taper)
# spacing between frequency points
df = 1. / (N * dt)
# grid of frequency values start at (N / 2) * df
f = (seq(N)  N / 2  1) * df
# radial frequency
omega = 2 * PI * f
# true spectrum (see example 1 on page 281)
spectrum = 1. / (1  2 * a*cos(omega) + a^2)
# only plot the positive frequencies
gxaxis("linear", 0., .5, 5, 5)
# log scale the y axis
gyaxis("log", 1e1, 1e2)
# title for the plot
gtitle("True and Estimated Spectrum")
# plot estimates with a plus sign
gplot(f, estimate, "plus")
# plot true spectrum with a solid line
gplot(f, spectrum, "solid")
gupdate
Output:

