knitr::opts_chunk$set( collapse = TRUE, comment = "#>" )
This vignette describes how to use the starmagarch packge to fit a STARMAGARCH($p,q,r,s$) model to data.
We start by simulating some data from a STARMAGARCH model. To do this, we must first specify the parameters and the neighbourhood system. The parameters should be a list containing the names mu
, phi
, theta
, omega
, alpha
, beta
, preferably in that order. The elements mu
and omega
are scalars (numeric), while the others should be matrices. The number of columns determines the maximum temporal lag and the number of rows the maximum spatial lag. These do not need to be equal for the differents parts of the model.
parameters <- list( mu = 5, phi = matrix(c(.8, .1), ncol = 1), #temporal order 1, spatial order 2 theta = matrix(c(-.42, -.2), ncol = 1), #temporal order 1, spatial order 2 omega = .9, alpha = matrix(c(.09, .13, .03, .07), ncol = 2), # temporal order 2, spatial order 2 beta = matrix(c(.5, .15), ncol = 1) #temporal order 1, spatial order 2 ) sum(parameters$alpha,parameters$beta)
The neighbourhood system is specified using the spdep
package. We have made a wrapper function that creates a three dimensional array where the first two dimensions give a neighbourhood matrix and the last is the order of the spatial lag. The first neighbourhood matrix (the lag 0 matrix) is an indentity matrix and all the remaining neighbourhood matrices are normalized so that each row (and column) add up to one. The spatial order of the neighbourhood system must be as large as the largest spatial order of the parameter list. In our example above, it must be at least 3 (lag 0 counts as 1). We also need to set a type
argument. This decides if a rook or queen neighbourhood should be enabled. The torus
argument decides whether or not the process should be circular, i.e. if the space should be "wrapped around a torus". The spatial sample size is set to $25\times 1$, i.e. a 1D space. The user may also specify other neighbourhood systems, but these are the automated ones.
# spatial dimension: m <- c(25, 1) library(starmagarch) W <- create.neighbourhood.array(m = m, sp = 3, type = "rook", torus =TRUE)
Now, we are ready to simulate a process. There are two calls, either you can specify the neighbourhood array in advance (as we just did) or the simulation function can create it for the specific simulation. The the spatial order is calculated from the parameters given. The reason for allowing these two options is that if you are to simulate multiple samples from the same process, it will be faster to not recreate the neighbourhood system every time. In this case, #1
is the faster, recommended option.
set.seed(1234) y <- simSTARMAGARCH(parameters, n = 2500, m = m, W=W, burnin = 500) # 1 #y <- simSTARMAGARCH(parameters, n = 2500, m = m, burnin = 500, # 2 #type= "rook", torus = TRUE)
library(ggplot2) library(reshape2) ggplot(data=melt(y), aes(Var2,Var1))+ geom_raster(aes(fill = value))+ scale_fill_gradient2(low="blue", mid="white", high = "red", midpoint = mean(y)) + xlab("Time") + ylab("Space") + theme_bw()
To Estimate a process, we must specify initial parameters. This also gives the order of the model you want to fit. Here we specify the model with correct order, but set some conservative values.
initial.parameters <- list( mu = mean(y), phi = matrix(c(.7, .01), ncol = 1), theta = matrix(c(.01, .01), ncol = 1), omega = 1, alpha = matrix(c(.01, .01, .01, .01), ncol = 2), beta = matrix(c(.01, .01), ncol = 1) )
The next step is to build a likelihood object from the TMB
package. We have created a template for STARMAGARCH and created a wrapper function that builds the model.
# Create likelihood object: # f <- CreateLikelihood(y, W=W, # init = apply(y,1,var), parameters=initial.parameters) # To save time, we use the true parameters as initial values: f <- CreateLikelihood(y, W=W, init = apply(y,1,var), parameters=parameters)
Here, f
is a list of functions containing the likelihood, gradient and hessian functions found by automatic differentiation. The init
argument allows the user to specify initial values for the ${\sigma_0^2(u)}$ process. This should be a non-negative vector of length prod(m)
. The other necessary initial values are set to zero.
Optimizing the likelihood is done using the nlminb
function. We have also here written a wrapper function to do this. This function also creates an object of the class starmagarch
, with classical generic functions attached to it, such as plot
, plot_garch
, summary
, AIC
, BIC
, residuals
, coef
, sigma
, fitted
, fittedgarch
, etc.
fit <- fitSTARMAGARCH(f, data = y, print = FALSE) summary(fit)
plot(fit) plot_garch(fit)
Comparing the results with the parameteres that generated the data:
round(rbind(unlist(parameters), round(coef(fit),4)), 4)
In some cases, we do not which to use a full scale ARMAGARCH model, but perhaps simply a GARCH model. This can be achieved by the following lines of code.
parameters <- list( mu = 0, phi = matrix(0, ncol = 1), #temporal order 1, spatial order 2 theta = matrix(0, ncol = 1), #temporal order 1, spatial order 2 omega = .9, alpha = matrix(c(.09, .13, .03, .07), ncol = 2), # temporal order 2, spatial order 2 beta = matrix(c(.5, .15), ncol = 1) #temporal order 1, spatial order 2 ) y <- simSTARMAGARCH(parameters, n = 2500, m = c(25,1), W = W) initial.parameters <- list( mu = mean(y), phi = matrix(c(.7, .01), ncol = 1), theta = matrix(c(.01, .01), ncol = 1), omega = 1, alpha = matrix(c(.01, .01, .01, .01), ncol = 2), beta = matrix(c(.01, .01), ncol = 1) ) f <- CreateLikelihood(y, W=W, init = apply(y,1,var), parameters = initial.parameters) fit <- fitSTARMAGARCH(f, data= y, print = FALSE) summary(fit) #> Estimates SD Zscore Pvalue #> mu 0.01272265 0.009833892 1.2937556 9.787493e-02 #> phi1 -0.37973220 0.282238217 -1.3454315 8.924296e-02 #> phi2 0.37536087 0.783248661 0.4792359 3.158854e-01 #> theta1 0.37367411 0.283705318 1.3171206 9.389911e-02 #> theta2 -0.36493089 0.791497331 -0.4610639 3.223764e-01 #> omega 0.98304477 0.066042385 14.8850586 2.060627e-50 #> alpha1 0.08262752 0.006244003 13.2331007 2.825359e-40 #> alpha2 0.14105291 0.011764199 11.9900142 2.004181e-33 #> alpha3 0.03092677 0.007764663 3.9830156 3.402315e-05 #> alpha4 0.06686734 0.014844508 4.5045170 3.326203e-06 #> beta1 0.51746041 0.033736850 15.3381362 2.126149e-53 #> beta2 0.10477532 0.048024420 2.1817093 1.456550e-02 #> #> #> Standardized residual standard error: 1 .
It seems that 'mu', 'phi1', 'phi2', 'theta1' and 'theta2' are not significantly different from zero. We can therefore fix these parameters to zero by using the 'map' argument.
# Use "map"" to fix parameters to zero: initial.parameters$mu <- 0 initial.parameters$phi <- matrix(0, ncol=1) initial.parameters$theta <- matrix(0, ncol=1) map <- list() map$mu <- as.factor(NA) map$phi <- as.factor(NA) map$theta <- as.factor(NA) dim(map$theta) <- dim(map$phi) <- c(1,1) map #> $mu #> [1] <NA> #> Levels: #> #> $phi #> [,1] #> [1,] <NA> #> Levels: #> #> $theta #> [,1] #> [1,] <NA> #> Levels: f <- CreateLikelihood(y, W=W, init = apply(y,1,var), parameters = initial.parameters, map = map) fit <- fitSTARMAGARCH(f, data= y, print = FALSE) summary(fit) #> Estimates SD Zscore Pvalue #> omega 0.98281266 0.066025733 14.885297 2.053285e-50 #> alpha1 0.08264733 0.006242327 13.239826 2.583416e-40 #> alpha2 0.14108928 0.011764874 11.992417 1.946865e-33 #> alpha3 0.03103404 0.007765063 3.996624 3.212613e-05 #> alpha4 0.06653875 0.014844846 4.482279 3.692497e-06 #> beta1 0.51727367 0.033725151 15.337920 2.133229e-53 #> beta2 0.10519927 0.048014155 2.191005 1.422570e-02 #> #> #> Standardized residual standard error: 1 . #> AIC: 148332.7 BIC: 148396 round(rbind(unlist(parameters)[-(1:3)], fit$coefficients),4) #> omega alpha1 alpha2 alpha3 alpha4 beta1 beta2 #> [1,] 0.9000 0.0900 0.1300 0.030 0.0700 0.5000 0.1500 #> [2,] 0.9828 0.0826 0.1411 0.031 0.0665 0.5173 0.1052
Equivalently, for a pure STARMA model:
parameters <- list( mu = 5, phi = matrix(c(.8, .1), ncol = 1), #temporal order 1, spatial order 2 theta = matrix(c(-.42, -.2), ncol = 1), #temporal order 1, spatial order 2 omega = 30, alpha = matrix(0, ncol = 1), # temporal order 2, spatial order 2 beta = matrix(0, ncol = 1) #temporal order 1, spatial order 2 ) y <- simSTARMAGARCH(parameters, n = 2500, m = c(25,1), W = W) initial.parameters <- list( mu = mean(y), phi = matrix(c(.7, .01), ncol = 1), theta = matrix(c(.01, .01), ncol = 1), omega = 1, alpha = matrix(c(.01, .01, .01, .01), ncol = 2), beta = matrix(c(.01, .01), ncol = 1) ) f <- CreateLikelihood(y, W=W, init = apply(y,1,var), parameters = initial.parameters) fit <- fitSTARMAGARCH(f, data= y, print = FALSE) summary(fit) #> Estimates SD Zscore Pvalue #> mu 5.0714676601 0.076651590 6.616259e+01 0.000000e+00 #> phi1 0.7977517884 0.005226702 1.526301e+02 0.000000e+00 #> phi2 0.1119027906 0.010694371 1.046371e+01 6.339603e-26 #> theta1 -0.4178454614 0.007787193 -5.365803e+01 0.000000e+00 #> theta2 -0.2233686304 0.016775178 -1.331543e+01 9.414894e-41 #> omega 29.8237638244 15.355764525 1.942187e+00 2.605724e-02 #> alpha1 0.0000000100 0.002948183 3.391920e-06 4.999986e-01 #> alpha2 0.0005736162 NaN NaN NaN #> alpha3 0.0001235004 0.002756660 4.480074e-02 4.821331e-01 #> alpha4 0.0001525837 0.009304228 1.639939e-02 4.934579e-01 #> beta1 0.0000000100 NaN NaN NaN #> beta2 0.0000000100 NaN NaN NaN #> #> #> Standardized residual standard error: 1 . #> AIC: 194676 BIC: 194784.5 # Use "map"" to fix parameters to zero: initial.parameters$alpha <- matrix(0, ncol=1) initial.parameters$beta <- matrix(0, ncol=1) map <- list() map$alpha <- as.factor(NA) map$beta <- as.factor(NA) dim(map$alpha) <- dim(map$beta) <- c(1,1) map #> $alpha #> [,1] #> [1,] <NA> #> Levels: #> #> $beta #> [,1] #> [1,] <NA> #> Levels: f <- CreateLikelihood(y, W=W, init = apply(y,1,var), parameters = initial.parameters, map = map) fit <- fitSTARMAGARCH(f, data= y, print = FALSE) summary(fit) #> Estimates SD Zscore Pvalue #> mu 5.0714315 0.076639614 66.17246 0.000000e+00 #> phi1 0.7977492 0.005225273 152.67130 0.000000e+00 #> phi2 0.1118951 0.010696428 10.46098 6.525081e-26 #> theta1 -0.4178449 0.007784903 -53.67374 0.000000e+00 #> theta2 -0.2233527 0.016775596 -13.31414 9.577871e-41 #> omega 29.8436922 0.168888961 176.70600 0.000000e+00 #> #> #> Standardized residual standard error: 1 . #> AIC: 194664 BIC: 194718.2 round(rbind(unlist(parameters)[1:6],coef(fit)), 4) #> mu phi1 phi2 theta1 theta2 omega #> [1,] 5.0000 0.8000 0.1000 -0.4200 -0.2000 30.0000 #> [2,] 5.0714 0.7977 0.1119 -0.4178 -0.2234 29.8437
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.