tfarima.R
In tfarima: Transfer Function and ARIMA Models

## ----setup, include=FALSE-------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 5,
  fig.height = 3.5,
  fig.align = "center",
  dpi = 120
)
library(tfarima)
data("BuildingMat")
data("seriesJ")
options(prompt = "R> ", continue = "+  ", width = 70, useFancyQuotes = FALSE)

## -------------------------------------------------------------------
# Creating an AR(1) process with mean 0.5 and innovation variance 1.5
ar1p <- um(ar = "1 - 0.8B", mu = 0.5, sig2 = 1.5)

# Displaying model properties
phi(ar1p)
roots(ar1p)
psi.weights(ar1p, 10)
autocorr(ar1p, 10)
spec(ar1p)[1:5, ]
display(ar1p)

## -------------------------------------------------------------------
set.seed(1234)
z <- sim(ar1p, n = 100)
ide(z)
ar1p.fit <- fit(ar1p, z)
diagchk(ar1p.fit)
p <- predict(ar1p.fit, n.ahead = 5)
plot(p, n.back = 20)

## -------------------------------------------------------------------
lp1 <- lagpol(param = c(theta = 0.5), p = 2)
lp2 <- lagpol(param = c(Theta1 = 1.2, Theta2 = -0.9), s = 12)
lp3 <- lagpol(param = c(theta = 0.8, Theta = 0.9),
              coef = c("theta", "Theta", "-theta*Theta"),
              lags = c(1, 12, 13))
lp3

## -------------------------------------------------------------------
i  <- lagpol0(list(1, c(1, 3), c(1, 12)), "i")
ma <- lagpol0("(1 - 0.8B)(1 - 0.7B3)(1 - 0.6B12)", "ma")
printLagpolList(i)
printLagpolList(ma)

## -------------------------------------------------------------------
ar1p <- um(ar = "(1 - 0.9B)")
ar1n <- um(ar = "(1 + 0.9B)")
ma1p <- um(ma = "(1 - 0.9B)")
ma1n <- um(ma = "(1 + 0.9B)")
ar2c <- um(ar = "(1 - 1.52B + 0.8B^2)")

display(list(ar1p, ar1n, ma1p, ma1n, ar2c),
        lag.max = 20, graphs = c("acf", "pacf", "spec"))

## -------------------------------------------------------------------
Z <- BuildingMat
ide(Z, transf = list(
  list(bc = TRUE),
  list(bc = TRUE, S = 1),
  list(bc = TRUE, d = 1, D = 1), 
  list(bc = TRUE, d = 1, D = 1, i = lagpol(s = 3, coef = "1"))
))

um1 <- um(Z, bc = TRUE,
          i  = list(1, c(1, 3), c(1, 12)),
          ma = list(theta = 1, Theta = c(1, 3), THETA = c(1, 12)))
diagchk(um1)

p <- predict(um1, n.ahead = 12, i = c(1, 12))
plot(p, n.back = 72)

## -------------------------------------------------------------------
uca1 <- as.ucarima(um1, i = 3)
uc1 <- wkfilter(uca1)
plot(uc1, main = "")

## -------------------------------------------------------------------
Y <- seriesJ$Y - mean(seriesJ$Y)
X <- seriesJ$X - mean(seriesJ$X)

umx <- um(X, ar = 3)
umy <- fit(umx, Y)
pccf(X, Y, um.x = umx, lag.max = 16)

tfx  <- tfest(Y, X, delay = 3, p = 2, q = 2, um.x = umx, um.y = umy)
tfmy <- tfm(Y, inputs = tfx, noise = um(ar = 2))
tfmy

## -------------------------------------------------------------------
Y <- window(BJsales, end = 140)
X <- window(BJsales.lead, end = 140)
umy <- um(Y, i = 1, ma = "1-0.620B", mu = 0.033, fit = FALSE)
umx <- um(X, i = 1, ma = "1 - 0.4483B", fit = FALSE)
tfx <- tf(X, delay = 3, w0 = 4.716, ar = "1-0.725B", um = umx)
tfmy <- tfm(output = Y, inputs = tfx, noise = umy, fit = FALSE)
plot(predict(tfmy, n.ahead = 10, ori = 89), n.back = 15, ylab = expression(Y[t]))

## -------------------------------------------------------------------
tfm1 <- outliers(um1, calendar = TRUE, easter = TRUE,
                 form = "wd", c = 4, p.value = 0.05, n.ahead = 24)
tfm1

## -------------------------------------------------------------------
trend <- um(i = 2, ma = c(theta = 1), sig2 = c(s2t = 0.1))
seas  <- um(i = "12", 
             ma = lagpol(coef = c("0.7944","0.6185","0.4699",
                                  "0.3465","0.2458","0.1658",
                                  "0.1041","0.0587","0.0275","0.0086")),
             sig2 = c(s2s = 0.01))
irreg <- um(sig2 = c(s2i = 1))

uca1 <- ucarima(Z, bc = TRUE,
                ucm = list(trend = trend, seas = seas, irreg = irreg))
rf1  <- as.um(uca1)
uca3 <- as.ucarima(rf1, i = 2, canonical = FALSE)
uca3