TOP:
SOLUTION:
The function to sample from an ARMA(p, q) process is as follows
my_arma = function(n, sigsq=1, mu=0, phi=numeric(0), theta=numeric(0), inits=NULL) { p = length(phi)
q = length(theta)
# Check that the process is stationary
if(p>0) {
roots = polyroot(c(1, -phi))
if(all(Mod(roots)>1)) message(“Process is stationary.”)
else message(“Process is not stationary.”)
}# Check that the process is invertible
if(q>0) {
roots = polyroot(c(1, theta))
if(all(Mod(roots)>1)) message(“Process is invertible.”)
else message(“Process is not invertible.”)
}
# Establish initial values
maxpq = max(c(p, q))
if(!is.null(inits)) {
if(length(inits)!=maxpq) stop(“If initial values are supplied,
there must be max(p, q) of them.”)
} else {
# Note: sets x[1], …, x[maxpq] = 0
inits = rep(0, maxpq)
}#
Simulate from model
x = numeric(n)
x[1:maxpq] = inits
eps = rnorm(n, 0, sqrt(sigsq))
for(t in (maxpq+1):n) {
x[t] = sum(phi * x[(t-1):(t-p)]) + sum(c(1, theta) * eps[t:(t-q)])
}r
eturn(x+mu)
}
We can sample a long realisation from this process as follows
(a) n = 1e4
sigsq = 0.2
mu = 1.5
phi = c(0.9, -0.5)
theta = -0.7
## Sample from process
x = my_arma(n, sigsq, mu, phi, theta)
## Process is stationary.
## Process is invertible.
The following function computes the autocovariance function
my_acf = function(x, lag.max, cor=FALSE) {
xbar = mean(x)
n = length(x)
gamma = numeric(lag.max+1)
for(k in 0:lag.max) { gamma[k+1] = sum((x[(1+k):n] – xbar) * (x[1:(n-k)] – xbar)) / n
}#
Allow the function to also compute the autocorrelation function
if(cor) gamma = gamma / gamma[1]
return(gamma)
}
(Note there is also a built-in R function – acf – that can be used). Now we can calculate
the requested moments as follows:
## Compute the sample mean
mean(x)
## [1] 1.496786
## Compute the sample autocovariance function for k=0,…,5
my_acf(x, 5)
## [1] 0.2680289835 0.0685127333 -0.0700936344 -0.0988780345 -0.0542393092
## [6] 0.0002972414
1. Consider the general ARIMA(p, d, q) process.
(a)Study the function from the TOP and then use it to give a new function which simulates from an ARIMA(p, d, q) process. This function should have d, the parameter vectors, and the length of the time series as arguments, and it should return the simulated time series. If you can, allow your function to handle any d. If not, write a function which handles the d = 0, d = 1 and d = 2 cases and returns with an error if d > 2(This second option will be marked for reduced credit.).
(b)
2. For an AR(p) model:
(a) Write a function which estimates the model parameters via the Yule-Walker equations.
(b) Test this function by simulating various AR processes and ensuring the original parameters are (approximately) recovered. Report the results of your tests. To simulate the
AR processes, you may like to use the function from the TOP.
3.
4. Consider the VARMA(p, q) process.
(a)Study the function from the TOP and then modify it to give a new function which simulates from a VARMA(p, q) process. This function should have the parameters, the length of the time series, and a set of initial values as arguments, and it should return the simulated time series. Tricky: make this function check if the entered parameters are such that the process is stationary and/or invertible.
(b)