R/pred.R

In PTSR: Positive Time Series Regression

#' @title Predict method for PTSR
#'
#' @description Predicted values based on ptsr object.
#'
#' @param object Object of class inheriting from \code{"ptsr"}
#' @param newdata A matrix with new values for the regressors.   If omitted
#' and \code{"xreg"} is present in the model, the fitted values are returned.
#' If the model does not include regressors, the functions will use
#' the value of \code{nnew}.
#' @param nnew number of out-of-sample forecasts required. If \code{newdata} is
#' provided, \code{nnew} is ignored.
#' @param ... further arguments passed to or from other methods.
#'
#' @details
#'   \code{predict.ptsr} produces predicted values, obtained by evaluating
#'   the regression function in the frame \code{newdata}.
#'
#'   If \code{newdata} is omitted the predictions are based on the data
#'   used for the fit.
#'
#'   For now, prediction intervals are not provided.
#'
#' @return A list with the following arguments
#'
#'  \item{series}{The original time series yt.}
#'
#'  \item{xreg}{The original regressors (if any).}
#'
#'  \item{fitted.values}{The in-sample forecast given by \eqn{\mu_t}.}
#'
#'  \item{etat}{In-sample values of \eqn{g(\mu[t])}.}
#'
#'  \item{error}{The error term}
#'
#'  \item{residuals}{The (in-sample) residuals, that is, the observed minus the predicted values.}
#'
#'  \item{forecast}{The predicted values for yt.}
#'
##' @export
predict.ptsr <-
  function(object, newdata, nnew = 0,...){

    out <- list()
    nms.out <- c("series", "xreg", "fitted.values", "residuals")

    if(missing(newdata)) newdata = NULL

    if(is.null(newdata) & nnew <= 0){
      #------------------------------------------------------
      # New data was not provided.
      # Extracting existing components and returning
      #------------------------------------------------------
      out[nms.out] <- object[nms.out]
    }else{

      if(is.null(newdata) & object$model$b > 0)
        stop("Please, provide the new values for the regressors")

      #------------------------------------------------------
      # New data was provided.
      # Making the necessary calculations
      #------------------------------------------------------
      xnew = NULL
      if(!is.null(newdata)){
        xnew = as.matrix(newdata)
        nnew = nrow(xnew)
      }

      temp <- .ptsr.predict(par = object$coefficients,
                            h = nnew, xreg = object$xreg,
														xnew = xnew, yt = object$series,
                            mut = object$fitted.values,
                            epst = object$residuals,
                            model = object$model)


      out[nms.out] <- object[nms.out]
      out$forecast <- temp
      out$xnew <- NULL
      if(object$model$b > 0) out$xnew <- xnew
    }
    out
  }


# Internal function.
# Calculates out-of-sample predicted values
.ptsr.predict <- function(par, h, xreg, xnew, yt, mut, epst, model) {

  n = length(yt)
  mut = c(mut, numeric(h))
  epst  = c(epst, numeric(h))

  p = model$p
  q = model$q
  a = model$a
  b = model$b
  arlag = model$arlag
  malag = model$malag
  g1 = model$g1
  g1.inv = model$g1.inv
  g2 = model$g2

  g2y = c(g2(yt), numeric(h))
  yt = c(yt, numeric(h))

  alpha = 0
  beta = 0
  l = 0
  if (a == 1) {alpha = par[1]; l = l + 1}
  if (b > 0) {beta = par[(l + 1):(l + b)]; l = l + b}
  if (p > 0) {
    if (is.null(arlag)) {
      phi = par[(l + 1):(l + p)]
      l = l + p
    }
    else {
      phi = numeric(p)
      phi[arlag] = par[(l + 1):(l + length(arlag))]
      l = l + length(arlag)
    }
  }
  if (q > 0) {
    if (is.null(arlag)) {
      theta = par[(l + 1):(l + q)]
      l = l + q
    }
    else{
      theta = numeric(q)
      theta[malag] = par[(l + 1):(l + length(malag))]
      l = l + length(malag)
    }
  }

  XB = numeric(n+h)
  if (!is.null(xnew)){
	  xtemp = rbind(as.matrix(xreg), as.matrix(xnew))
		XB = xtemp %*% beta
	}

  for (t in 1:h) {
    ls = alpha + XB[n+t]
    if (p > 0) {
      xr = 0
      if (model$xregar) xr = XB[(n + t) - c(1:p)]
      temp = sum(phi * (g2y[(n + t) - c(1:p)] - xr))
      ls = ls + temp
    }
    if (q > 0) ls = ls + sum(theta * epst[(n + t) - (1:q)])
    mut[n + t] = g1.inv(ls)
    yt[n + t] = mut[n + t]
    g2y[n+t] = g2(mut[n + t])
  }

  ynew = yt[(n + 1):(n + h)]

  invisible(ynew)
}