Nothing
####################################################################
# This program estimates the parameters of an GLARMA(p,q) model
# for the Poisson Distribution by using Fisher scoring or
# Newton-Raphson iterations.
#
# Score type (GAS) Residuals are used in the state equation.
#
# User Manual:
# Dunsmuir (2010) "R software for ?fitting observations driven
# regression models for univariate time series"
#
# Contributers: See Dunsmuir(2010).
#
# Last Modified: March 29, 2011.
#
####################################################################
glarmaPoissonScore <- function(y, X, offset = NULL, delta, phiLags, thetaLags,
method = "FS"){
r <- ncol(X)
n <- length(y)
p <- length(phiLags)
q <- length(thetaLags)
beta <- delta[1:r]
phi <- delta[(r + 1):(r + p)]
theta <- delta[(r + p + 1):(r + p + q)]
mpq <- 0
if ((p + q) > 0) {
mpq <- max(phiLags[p], thetaLags[q])
}
nmpq <- n + mpq
s <- r + p + q
e <- array(0, nmpq)
Z <- array(0, nmpq)
W <- array(0, nmpq)
mu <- array(0, nmpq)
e.d <- array(0, c(s, nmpq))
Z.d <- array(0, c(s, nmpq))
W.d <- array(0, c(s, nmpq))
if (method == "NR") {
e.dd <- array(0, c(s, s, nmpq))
Z.dd <- array(0, c(s, s, nmpq))
W.dd <- array(0, c(s, s, nmpq))
}
if(is.null(offset)) eta<-X %*% beta else eta<- X %*% beta + offset
ll <- 0
ll.d <- matrix(0, ncol = 1, nrow = s)
ll.dd <- matrix(0, ncol = s, nrow = s)
## Start loop on time
for (time in 1:n) {
tmpq <- time + mpq
if (p > 0) {
Z.d[(r + 1):(r + p), tmpq] <- Z[tmpq - phiLags] + e[tmpq - phiLags]
if (method == "NR") {
Z.dd[(r + 1):(r + p), , tmpq] <- t((Z.d + e.d)[, (tmpq - phiLags)])
Z.dd[, (r + 1):(r + p), tmpq] <- Z.dd[, (r + 1):(r + p), tmpq] +
(Z.d + e.d)[, (tmpq - phiLags)]
}
for (i in 1:p) {
Z[tmpq] <- Z[tmpq] + phi[i] * (Z + e)[tmpq - phiLags[i]]
Z.d[, tmpq] <- Z.d[, tmpq] + phi[i] * (Z.d[, tmpq - phiLags[i]] +
e.d[, tmpq - phiLags[i]])
if (method == "NR") {
Z.dd[, , tmpq] <- Z.dd[, , tmpq] + phi[i] *
(Z.dd[, , tmpq - phiLags[i]] +
e.dd[, , tmpq - phiLags[i]])
}
}
}
if (q > 0) {
Z.d[(r + p + 1):(r + p + q), tmpq] <- e[tmpq - thetaLags]
if (method == "NR") {
Z.dd[(r + p + 1):(r + p + q), , tmpq] <- Z.dd[(r + p + 1):
(r + p + q), , tmpq] +
t(e.d[, tmpq - thetaLags])
Z.dd[, (r + p + 1):(r + p + q), tmpq] <- Z.dd[, (r + p + 1):
(r + p + q), tmpq] +
e.d[, tmpq - thetaLags]
}
for (i in 1:q) {
Z[tmpq] <- Z[tmpq] + theta[i] * e[tmpq - thetaLags[i]]
Z.d[, tmpq] <- Z.d[, tmpq] + theta[i] * e.d[, tmpq - thetaLags[i]]
if (method == "NR") {
Z.dd[, , tmpq] <- Z.dd[, , tmpq] + theta[i] *
e.dd[, , tmpq - thetaLags[i]]
}
}
}
W[tmpq] <- eta[time] + Z[tmpq]
W.d[, tmpq] <- matrix(c(X[time, ], rep(0, p + q)), ncol = 1) +
Z.d[, tmpq]
if (method == "NR") {
W.dd[, , tmpq] <- Z.dd[, , tmpq]
}
## Poisson mean
mu[tmpq] <- exp(W[tmpq])
## Score Residuals
e[tmpq] <- (y[time] - mu[tmpq])/mu[tmpq]
## update derivates of Score Residuals for Poisson
e.d[, tmpq] <- -(1 + e[tmpq]) * W.d[, tmpq]
if (method == "NR") {
e.dd[, , tmpq] <- -(1 + e[tmpq]) * W.dd[, , tmpq] - (e.d[, tmpq] *
W.d[, tmpq]) %o% W.d[, tmpq]
}
# Note the above could be simplified to:
# e.d[, tmpq] <- -(y[time]/mu[tmpq]) * W.d[, tmpq]
# if (method == "NR") {
# e.dd[, , tmpq] <- -(y[time]/mu[tmpq]) * ( W.dd[, , tmpq] -
# W.d[, tmpq]) %o% W.d[, tmpq]
# }
## update likelihood and derivatives for Poisson
ll <- ll + y[time] * W[tmpq] - mu[tmpq] - log(factorial(y[time]))
ll.d <- ll.d + (y[time] - mu[tmpq]) * W.d[, tmpq]
if (method == "FS") {
ll.dd <- ll.dd - mu[tmpq] * W.d[, tmpq] %o% W.d[, tmpq]
}
if (method == "NR") {
ll.dd <- ll.dd + (y[time] - mu[tmpq]) * W.dd[, , tmpq] - mu[tmpq] *
W.d[, tmpq] %o% W.d[, tmpq]
}
}
## end loop on time.
list(delta = delta, ll = ll, ll.d = ll.d, ll.dd = ll.dd, eta = eta,
W = W[mpq + 1:n], e = e[mpq + 1:n], mu = mu[mpq + 1:n],
fitted.values = mu[mpq + 1:n])
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.