Nothing
##' @title arma.forecast
##'
##' @description Forecasting of (multivariate) time series of
##' using marima type model.
##'
##' @param series matrix holding the kvar-variate timeseries.
##' The series is assumed to have the same format
##' as the timeseries analysed by marima BEFORE differencing (if
##' differencing was used via define.dif)
##' (the length, though, does not need to be the same but can be shorter
##' or longer). Results
##' from estimating the model (for the differenced data, if used) are assumed
##' to be saved in the input-object 'marima' (see 'usage') by marima.
##'
##' The series is assumed to have the total length=(nstart+nstep) (but it
##' may be longer. In any case the forecasting is starting from nstart
##' continuing to nstart+nstep. Future values already present or initialised,
##' for example, as NAs are overwritten with the forecasted values.)
##'
##' An example of a series prepared for forcasting is in the marima library:
##' 'data(austr)': (see below, the example).
##'
##' If future (independent) x-values for the forecasting are to be used
##' these values must be supplied in 'series' at the proper places before
##' calling 'arma.forecast(...)' (that is except the x-value(s)
##' corresponding to the last prediction).
##'
##' @param marima the object holding the marima results to be used for
##' the forecasting, that is an output object created by marima.
##'
##' If the ar- and/or the ma-model do not include a leading unity matrix
##' this is automatically taken care of in the function (in that case the
##' dimensions of the model arrays used will be, respectively,
##' (kvar, kvar, p+1) and (kvar, kvar, q+1)) after inserting the leading
##' unity matrix (if the object 'marima' was produced by marima, this
##' will automatically be OK.
##'
##' @param nstart starting point for forecasting (1st forecast values
##' will be for time point t = nstart+1).
##'
##' @param nstep length of forecast (forecasts will be for time points
##' nstart+1,...,nstart+nstep).
##'
##' @param dif.poly (most often) output from the function define.dif holding
##' the ar-representation of the differencing polynomial
##' (define.dif$dif.poly).
##' If a differenced timeseries was analysed by marima
##' the forecast-variance/covariance matrices are calculated for the
##' aggregated (original) timeseries if 'dif.poly' is specified. If not,
##' the forecast-variance/covariance matrices are calculated for the
##' differenced time series. If forecasting is wanted for the original
##' (not differenced) time series the 'dif.poly' created by define.dif
##' must be specified.
##'
##' @param check If check=TRUE (default) various checks and
##' printouts are carried out.
##'
##' @return forecasts = forecasted values following the nstart first values
##' of the input series (at time points 'nstart+1,...,nstart+nstep').
##' The forecasted values will be (over-) written in the input series at
##' the proper future positions (if relevant).
##'
##' @return residuals = corresponding residuals for input series followed by
##' nstep future residuals (all=0).
##'
##' @return prediction.variances = (kvar, kvar, nstep) array containing
##' prediction covariance matrices corresponding to the nstep forecasts.
##' @return nstart = starting point for prediction (1st prediction at point
##' nstart+1).
##' @return nstep = length of forecast
##'
##' @examples
##'
##' library(marima)
##' data(austr)
##' series<-austr
##' Model5 <- define.model(kvar=7, ar=1, ma=1, rem.var=1, reg.var=6:7)
##' Marima5 <- marima(ts(series[1:90, ]), Model5$ar.pattern, Model5$ma.pattern,
##' penalty=1)
##'
##' nstart <- 90
##' nstep <- 10
##' cat("Calling arma.forecast.\n")
##' cat("In the example the input series is dim(length,kvar).\n")
##' cat("and of type ts() (timeseries) for illustration. \n")
##' Forecasts <- arma.forecast(series=ts(series), marima=Marima5,
##' nstart=nstart, nstep=nstep )
##' Year<-series[91:100,1]
##' One.step <- Forecasts$forecasts[, (nstart+1)]
##' One.step
##' Predict <- Forecasts$forecasts[ 2, 91:100]
##' Predict
##' stdv<-sqrt(Forecasts$pred.var[2, 2, ])
##' upper.lim=Predict+stdv*1.645
##' lower.lim=Predict-stdv*1.645
##' Out<-rbind(Year, Predict, upper.lim, lower.lim)
##' print(Out)
##' # plot results:
##' plot(series[1:100, 1], Forecasts$forecasts[2, ], type='l', xlab='Year',
##' ylab='Rate of armed suicides', main='Prediction of suicides by firearms',
##' ylim=c(0.0, 4.1))
##' lines(series[1:90, 1], series[1:90, 2], type='p')
##' grid(lty=2, lwd=1, col='black')
##' Years<-2005:2014
##' lines(Years, Predict, type='l')
##' lines(Years, upper.lim, type='l')
##' lines(Years, lower.lim, type='l')
##' lines(c(2004.5, 2004.5), c(0.0, 2.0), lty = 2)
##'
##' @importFrom graphics grid plot
##'
##' @export
arma.forecast <- function( series=NULL, marima=NULL, nstart=NULL,
nstep=1, dif.poly=NULL , check=TRUE ) {
if(check) {
if( is.ts(series) ){
series <- as.matrix(series)
cat("\n Input series is type ts (timeseries).
It will be changed to matrix(...) using as.matrix(...). \n" )
}
}
Y <- series
d <- dim(Y)
if (d[1] > d[2]) {
if(check){ cat("\n", "Data matrix is to be transposed. Done! \n") }
Y <- t(Y)
d <- dim(Y)
}
cat(d[1], " Variables in series", ", with required length = ",
(nstart+nstep),". \n")
Y <- Y[ , 1:(nstart + nstep) ]
if (is.null(dif.poly)) {
dif.poly <- array(diag(d[1]), dim = c( d[1], d[1], 1 ))
}
dif.poly <- check.one( dif.poly )
colnames(Y) <- c(1:d[2])
means <- marima$mean.pattern
averages <- marima$averages
Constant <- marima$Constant
if (is.null(nstart)) {
nstart <- d[2]
}
Names <- rownames(marima$y.analysed)
rownames(Y) <- Names
AR <- pol.mul( marima$ar.estimate, dif.poly,
L <- ( dim(marima$ar.estimate)[3] + dim(dif.poly)[3]-2 ) )
MA <- marima$ma.estimate
p <- dim(AR)[3]
q <- dim(MA)[3]
if(check) {
cat("\n print AR-model: \n")
print(AR)
cat("\n print MA-model: \n")
print(MA)
}
Rand <- rep(TRUE, d[1])
Regr <- !Rand
sig2 <- marima$resid.cov
# print(sig2,digits=2)
for (i in 1:d[1]) {
if (sum(abs(AR[i, , ])) + sum(abs(MA[i, , ])) == 2) {
Rand[i] <- FALSE
sig2[i, ] <- 0
sig2[, i] <- 0
}
if (sum(abs(AR[, i, ])) + sum(abs(MA[, i, ])) > 2) {
Regr[i] <- TRUE
}
}
L <- nstart + nstep
if (L > d[2]) {
cat("Input series length =", d[2], " but (nstart+nstep)=", L, "\n")
}
for (i in 1:d[1]) {
# cat(Rand[i],Regr[i],'\n')
if (!Rand[i] & Regr[i]) {
cat("variable no.", i, "not random and regressor.\n")
zy <- Y[i, (nstart:(L-1))]
# zy[5]<-NaN
if (is.na(sum(zy)) | is.nan(sum(zy))) {
cat("Variable no.", i, ": illegal (NaN or NA) in future values. \n")
}
if (!is.na(sum(zy)) & !is.nan(sum(zy))) {
cat("Variable no.", i, ": future values seem OK.\n")
}
}
}
# cat('Calling series : \n') print(Y[,(d[2]-19):d[2]])
forecast <- arma.forecastingY(series = Y, ar.poly = AR,
ma.poly = MA, nstart = nstart, nstep = nstep, Rand = Rand,
Regr = Regr, marima = marima)
pred.var <- forec.var(marima, nstep = nstep, dif.poly )$pred.var
return(list(nstart = nstart, nstep = nstep,
forecasts = forecast$estimates,
residuals = forecast$residuals,
pred.var = pred.var, averages = averages,
mean.pattern = means))
}
arma.forecastingY <- function(series = NULL, ar.poly = NULL,
ma.poly = NULL, nstart = NULL, nstep = NULL, Rand = NULL,
Regr = NULL, marima = NULL) {
"[" <- function(x, ...) .Primitive("[")(x, ..., drop = FALSE)
series <- as.matrix(series)
d <- dim(series)
if (d[1] > d[2]) {
series <- t(series)
}
d <- dim(series)
k <- d[1]
N <- d[2]
kvar <- k
means <- marima$mean.pattern
averages <- marima$averages
Constant <- marima$Constant
cat("Constant =",marima$Constant,"\n")
# cat("dimensions of time series", k, N, "\n")
ar <- dim(ar.poly)
ma <- dim(ma.poly)
yseries <- matrix(c(series), nrow = d[1])
yseries[is.na(yseries) | is.nan(yseries)] <- 0
residuals <- yseries * 0
forecasts <- residuals
kvar <- ar[1]
extra <- ar[3]
su0 <- matrix(0, nrow=kvar, ncol=1)
if( extra > 1 ){
for (i in 1:extra){
forecasts[, i] <- yseries[, i]
residuals[, i] <- yseries[, i] - forecasts[, i]
}
}
for (i in extra:N){
suar <- su0
if(ar[3] > 1) {
for (j in 2:ar[3]) {
suar <- suar + matrix(ar.poly[, , j], nrow = kvar) %*%
matrix(yseries[, (i+1-j)], ncol=1)
if (i < 0)
cat("i, suar= ", i, "\n")
}
}
suma <- su0
if (ma[3] > 1) {
for (j in 2:ma[3]) {
if( i+1-j > 0 ) {
suma <- suma + matrix(ma.poly[, , j], nrow = kvar) %*%
matrix(residuals[, (i + 1 - j)], ncol = 1)
}
}
}
# cat("suar = ",suar,"\n"); print(suar)
# cat("suma = ",suma,"\n"); print(suma)
# cat("Constant = ",Constant,"\n");
# cat("Constant dimension=", dim(Constant),"\n"); print(Constant)
est <- suma - suar + Constant
forecasts[, i] <- est
# cat('ar- & ma-polynomier \n')
# print(ar.poly[,,2])
# print(ma.poly[,,2])
if (i <= nstart) {
for (j in 1:k) {
if (Rand[j]) {
residuals[j, i] <- yseries[j, i] - est[j]
}
}
}
if (i > nstart) {
for (j in 1:k) {
if (Rand[j]) {
yseries[j, i] <- est[j]
residuals[j, i] <- 0
}
}
}
} # End i-loop here
return(list(estimates = forecasts, residuals = residuals))
}
##' @title forec.var
##'
##' @description Function for calculation of variances of nstep forecasts
##' using a marima type model.
##'
##' @param marima marima object (cov.u and ar.estimates and
##' ma.estimates are used)
##' @param nstep length of forecast
##'
##' @param dif.poly autoregressive representation of differencing
##' polynomial as constructed by the function
##' define.dif(...) when the time series is differenced
##' (if so) before being analysed by marima.
##'
##' @return pred.var = variance-covariances for nstep forecasts
##' (an array with dimension (kvar, kvar, nstep).
##'
##' @return rand.shock = corresponding random shock representation
##' of the model used.
##'
forec.var <- function(marima, nstep = 1, dif.poly = NULL) {
cat("Calculation of forecasting variances. \n")
sig2 <- marima$resid.cov
ar.poly <- marima$ar.estimates
ma.poly <- marima$ma.estimates
d <- dim(sig2)
if(is.null(dif.poly)) {
dif.poly <- diag(d[1])
}
dif.poly <- check.one(dif.poly)
kvar = dim(sig2)[1]
ar.poly<-pol.mul(ar.poly, dif.poly, L=nstep+1)
if (is.null(ma.poly)) {
ma.poly <- diag(kvar)
}
ar.poly <- check.one(ar.poly)
ma.poly <- check.one(ma.poly)
if (nstep < 1) {
nstep <- 1
}
xsi.poly <- rand.shock(ar.poly, ma.poly, nstep)
var <- xsi.poly
# cat('nsteps=',1:nstep,'\n')
for (i in 1:nstep) {
# cat('i=',i,'\n')
var[, , i] <- matrix(xsi.poly[, , i], nrow = kvar) %*%
sig2 %*% t(matrix(xsi.poly[, , i], nrow = kvar))
if (i > 1) {
var[, , i] <- var[, , i] + var[, , (i - 1)]
}
}
d <- dim(var)[3] - 1
# cat('d=',d,'\n')
# cat('xsi.poly=',1:nstep,'\n')
# print(round(xsi.poly,4))
# cat('var=','\n') print(round(var[,,1:d],4))
return(list(pred.var = var[, , 1:d], rand.shock = xsi.poly[, , 1:d]))
}
##' @title arma.filter
##'
##' @description Filtering of (kvar-variate) time series with marima
##' type model.
##'
##' Calculation of residuals and filtered values of timeseries using
##' a marima model.
##'
##' @param series matrix holding the kvar by n multivariate timeseries
##' (if (kvar > n) the series is transposed and a warning is given).
##'
##' @param ar.poly (kvar, kvar, p+1) array containing autoregressive matrix
##' polynomial model part. If the filtering is to be performed for
##' undifferenced data when the analysis (in marima) was done for differenced
##' data, the input array ar.poly should incorporate the ar-representation
##' of the differensing operation (using, for example:
##' ar.poly <- pol.mul(ar.estimate, dif.poly,
##' L = ( dim(ar.estimates)[3]+dim(dif.poly)[3])), where 'dif.poly'
##' was obtained when differencing the time series (using define.dif)
##' before analysing it with marima (giving the ar.estimate) .
##'
##' @param ma.poly (kvar, kvar, q+1) array containing moving average matrix
##' polynomial model part.
##'
##' If a leading unity matrix is not included in the ar- and/or the ma-part
##' of the model this is automatically taken care of in the function
##' (in that case the dimensions of the model arrays used in arma.filter()
##' are, respectively, (kvar, kvar, p+1) and (kvar, kvar, q+1)).
##'
##' @param means vector (length = kvar) indicating whether means are
##' subtracted or not (0/1). Default : means = 1 saying that all means
##' are subtracted (equivalent to means = c(1, 1, ..., 1)).
##'
##' @return estimates = estimated values for input series
##'
##' @return residuals = corresponding residuals. It is noted that the
##' residuals computed by arma.filter may deviate slightly from the
##' marima-residuals (which are taken from the last repeated regression
##' step performed). The residuals computed by arma.filter are constructed
##' by filtering (successive use of the arma model) and using a heuristic
##' method for the first residuals.
##'
##' @return averages = averages of variables in input series
##'
##' @return mean.pattern = pattern of means as used in filtering
##'
##' @examples
##'
##' library(marima)
##' data(austr)
##' series<-t(austr)[,1:90]
##' # Define marima model
##' Model5 <- define.model(kvar=7,ar=1,ma=1,rem.var=1,reg.var=6:7)
##'
##' # Estimate marima model
##' Marima5 <- marima(series,Model5$ar.pattern,Model5$ma.pattern,penalty=1)
##'
##' # Calculate residuals by filtering
##' Resid <- arma.filter(series, Marima5$ar.estimates,
##' Marima5$ma.estimates)
##'
##' # Compare residuals
##'
##' plot(Marima5$residuals[2, 5:90], Resid$residuals[2, 5:90],
##' xlab='marima residuals', ylab='arma.filter residuals')
##'
##' @export
##'
arma.filter <- function(series = NULL, ar.poly = array(diag(kvar),
dim = c(kvar, kvar, 1)), ma.poly = array(diag(kvar),
dim = c(kvar, kvar, 1)), means = 1) {
"[" <- function(x, ...) .Primitive("[")(x, ..., drop = FALSE)
cat("arma.filter is being called \n")
if (is.null(series)) {
stop("No input series specified in input to arma.filter \n")
}
series <- as.matrix(series)
d <- dim(series)
if (d[1] > d[2]) {
cat("Warning: Input series is transposed to be a k x n series. \n")
cat("Output (estimated values and residuals) \n")
cat("will be organised the same way (k x n). \n")
series <- t(series)
}
vmeans <- means
# cat("vmeans,means=", vmeans, means, "\n")
d <- dim(series)
kvar <- d[1]
averages <- rep(0, kvar)
means <- rep(1, kvar)
if (length(vmeans) == 1) {
if (vmeans != 1) {
means <- rep(0, kvar)
}
}
for (i in 1:kvar) {
averages[i] <- mean(series[i, ])
if (means[i] == 1) {
series[i, ] <- series[i, ] - averages[i]
}
}
cat("indicators for means=", means, "\n")
# cat(averages, "\n")
# ar.poly <- check.one(ar.poly)
# ma.poly <- check.one(ma.poly)
m <- dim(ar.poly)
ma <- dim(ma.poly)
# cat("m = ",m," \n")
# cat("ma = ",ma,"\n")
su0 <- matrix(0, nrow = kvar, ncol = 1)
extra <- 2 * max(m[3], ma[3])
extray <- matrix(0, nrow = kvar, ncol = extra)
for (i in 1:kvar) {
if (means[i] != 1) {
extray[i, ] <- extray[i, ] + averages[i]
}
}
# cat('extra=',extra,'\n') print(extray)
yseries <- cbind(extray, series)
estimates <- yseries * 0
residuals <- estimates
# cat('Series= \n') print(averages) print(round(yseries[,1:20],2))
# print(round(short.form(ar.poly,leading=F),2))
s <- dim(yseries)
cat(" dim(yseries)", s, "\n")
for (i in (extra/2 + 1):s[2]) {
sur <- su0
if (m[3] > 1) {
for (j in 2:m[3]) {
# cat('sur=',sur,'\n') cat('m,i,j=',m,i,j,'\n')
# cat('dimensions=',dim(ar.poly[,,j]),dim(yseries[,(i+1-j)]),
# '\n')
# print(ar.poly[,,j]) print(yseries[,(i+1-j)])
sur <- sur + matrix(ar.poly[, , j], nrow = kvar) %*%
matrix(yseries[, (i + 1 - j)], nrow = kvar)
}
}
# cat(i,j,(i+1-j),sur,'\n') cat('y=',yseries,'\n')
suma <- su0
if (ma[3] > 1) {
for (j in 2:ma[3]) {
suma <- suma + matrix(ma.poly[, , j], nrow = kvar) %*%
matrix(residuals[, (i + 1 - j)], nrow = kvar)
}
}
# cat(i,j,(i+1-j),suma,'\n')
su <- -sur + suma
estimates[, i] <- su
residuals[, i] <- yseries[, i] - estimates[, i]
}
estimates <- estimates[, (extra + 1):(extra + d[2])]
residuals <- residuals[, (extra + 1):(extra + d[2])]
for (i in 1:kvar) {
if (means[i] == 1) {
estimates[i, ] <- estimates[i, ] + averages[i]
}
}
# cat('dimensions',dim(yseries),dim(series),dim(residuals),'\n')
return(list(estimates = estimates,
residuals = residuals,
averages = averages,
mean.pattern = means))
}
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.