library("UComp")
library("ggplot2")
library("gridExtra")
library("forecast")
library("autostsm")
library("zoo")
library("tsoutliers")
library("rucm")
library("data.table")
theme_set(theme_minimal())

##################################
# Auxiliar functions:
# rowMedians: median of row of matrix
# ForecFun: runs all models in paper for one forecasting origin
# metrics:  computes error metrics for one model and one forecasting origin
# rolling:  runs a rolling window experiment for many forecasting origins
##################################
rowMedians <- function(x, na.rm = FALSE, ...) {
    apply(x, 1, median, na.rm = na.rm, ...)
}
forecFun <- function(y) {
    # Function to forecast all models
    n <- length(y)
    y <- log(y)
    # ETS
    mETS <- ets(y)
    fETS <- forecast(mETS, h = 12)$mean
    # ARIMA
    mARIMA <- auto.arima(y)
    fARIMA <- forecast(mARIMA, h = 12)$mean
    # UC
    mUC <- UCmodel(y, h = 12)
    fUC <- mUC$yFor
    # BSM
    mBSM <- UCmodel(y, h = 12, model = "llt/equal/arma(0,0)")
    fBSM <- mBSM$yFor
    # Seasonal naive
    fNaive <- y[(n - 11):n]
    fNaive <- ts(fNaive, start = start(fUC), frequency = frequency(y))
    # auto_stsm
    x <- data.table(data.frame(date = as.Date(y), y = as.matrix(y)))
    stsm <- stsm_estimate(x, freq = 12, cycle = FALSE)
    stsm_fc <- stsm_forecast(stsm, y = x, n.ahead = 12)
    fAuto <- ts(stsm_fc$fitted[(n + 1):(n + 12)], start = start(fUC), frequency = frequency(y))
    # Mean
    fMean <- rowMeans(cbind(fETS, fARIMA, fUC, fBSM, fNaive))
    # Median
    fMedian <- rowMedians(cbind(fETS, fARIMA, fUC, fBSM, fNaive))
    return(exp(cbind(fETS, fARIMA, fUC, fBSM, fAuto, fMean, fMedian, fNaive)))
}
metrics <- function(insample, outsample, forecasts) {
    # Compute sMAPE and MASE for a single forecasting origin, method and horizon
    absE <- abs(outsample - forecasts)
    sMAPE <- mean(absE * 200 / (abs(outsample) + abs(forecasts)), na.rm = TRUE)
    n <- length(insample)
    MASEDenominator <- abs(insample[(13):n] - insample[1:(n - 12)])
    Denominator <- mean(MASEDenominator, na.rm = TRUE)
    MASE <- mean(absE, na.rm = TRUE) / Denominator
    return(round(cbind(sMAPE, MASE), 3))
}
rolling <- function(y, orig) {
    # Function to compute error metrics in a rolling experiment
    n <- length(y)
    origs <- orig:(n - 12)
    nOr <- length(origs)
    methods <- c("ETS", "ARIMA", "UC", "BSM", "Auto", "Mean", "Median", "Naive")
    nMethods <- length(methods)
    sMAPEaux <- array(NA, c(12, nOr, nMethods))
    MASEaux <- sMAPEaux
    # Loop for forecasting origin
    for (i in 1:nOr) {
        cat(paste(i, "of", nOr, "\n"))
        yi <- subset(y, end = origs[i])
        outi <- forecFun(yi)
        # Loop for horizon
        for (j in 1:12) {
            # Loop for method
            for (k in 1:nMethods) {
                aux <- metrics(yi, subset(y, start = origs[i] + 1, end = origs[i] + j), outi[1:j, k])
                sMAPEaux[j, i, k] <- aux[1]
                MASEaux[j, i, k] <- aux[2]
            }
        }
        # print(MASEaux[, i, ])
    }
    sMAPE <- matrix(NA, 12, nMethods)
    MASE <- sMAPE
    for (i in 1:nMethods) {
        sMAPE[, i] <- rowMeans(sMAPEaux[,,i])
        MASE[, i] <- rowMeans(MASEaux[,,i])
    }
    colnames(sMAPE) <- methods
    colnames(MASE) <- methods
    return(list(sMAPE = sMAPE, MASE = MASE))
}

#####################################
# Example 1: Air Passengers data from Box and Jenkins (2015)
#####################################
y <- log(AirPassengers)
autoplot(AirPassengers, ylab = "Thousands")
# Estimating Basic Structural Model
m <- UC(y, model = "llt/equal/arma(0,0)")
summary(m)
# Full algorithm
m <- UC(y)
# Stepwise with unit roots
m <- UC(y, stepwise = TRUE, tTest = TRUE)
# Stepwise
m <- UC(y, stepwise = TRUE)
m
# Plotting components
t <- time(y)
# comp <- ts(m$comp[1:length(t), ], start = start(y), frequency = 12)
comp <- subset(m$comp, end = 144)
comp[is.na(comp)] <- 0
p1 <- ggplot() + 
    geom_line(data = y, aes(x = t, y = y)) +
    geom_line(data = comp[, 1], aes(x = t, y = comp[,1]), color = "red") +
    ylab("Trend") +
    theme(axis.title.x=element_blank(),
          axis.text.x=element_blank(),
          axis.ticks.x=element_blank())
p2 <- ggplot() + 
    geom_line(data = comp[, 3], aes(x = t, y = comp[, 3])) +
    xlab("") + ylab("Seasonal")+
    theme(axis.title.x=element_blank(),
          axis.text.x=element_blank(),
          axis.ticks.x=element_blank())
p3 <- ggplot(data = comp[, 4], aes(x = t, y = comp[, 4])) + 
    geom_bar(stat = "identity", width = 0.01, fill = "#000000") + 
    xlab("Years") + ylab("Irregular")
grid.arrange(grobs = list(p1, p2, p3),
    widths = 1,
    layout_matrix = rbind(1, 1, 2, 3))
# dev.print(device = pdf, file = "airpas.pdf")
# Plotting forecasts
m <- UC(subset(y, end = 132), h = 12)
py <- subset(y, start = 120)
forec <- predict(m)
ggplot() + 
  geom_line(data = py, aes(x = time(py), y = py), linetype = "solid") +
  geom_line(data = forec[, 1], aes(x = time(forec), y = forec[, 1]), linetype = "dashed") +
  geom_line(data = forec[, 2], aes(x = time(forec), y = forec[, 2]), linetype = "dotted", color = "red") +
  geom_line(data = forec[, 3], aes(x = time(forec), y = forec[, 3]), linetype = "dotted", color = "red") + 
  geom_vline(xintercept = 1960) +
  xlab("Years") + ylab("Air passengers (log scale)")
# dev.print(device = pdf, file = "airpasFore.pdf")
# Table 5 in paper
firstRow <-     UC(y)$criteria
secondRow <-    UC(y, criterion = "bic")$criteria
secondRowBis <- UC(y, criterion = "aicc")$criteria
thirdRow <-     UC(y, model = "llt/equal/arma(0,0)")$criteria
print(rbind(firstRow, secondRow, secondRowBis, thirdRow))
# Table 5 in paper stepwise
firstRow <-     UC(y, stepwise = TRUE)$criteria
secondRow <-    UC(y, criterion = "bic", stepwise = TRUE)$criteria
secondRowBis <- UC(y, criterion = "aicc", stepwise = TRUE)$criteria
thirdRow <-     UC(y, model = "llt/equal/arma(0,0)", stepwise = TRUE)$criteria
print(rbind(firstRow, secondRow, secondRowBis, thirdRow))
# Table 5 in paper stepwise with unit roots tests
firstRow <-     UC(y, stepwise = TRUE, tTest = TRUE)$criteria
secondRow <-    UC(y, criterion = "bic", stepwise = TRUE, tTest = TRUE)$criteria
secondRowBis <- UC(y, criterion = "aicc", stepwise = TRUE, tTest = TRUE)$criteria
thirdRow <-     UC(y, model = "llt/equal/arma(0,0)", stepwise = TRUE, tTest = TRUE)$criteria
print(rbind(firstRow, secondRow, secondRowBis, thirdRow))
# Rolling forecasting exercise from observation 108
errors <- rolling(AirPassengers, 108)
print(errors$sMAPE)
print(errors$MASE)

#####################################
# Example 2: Sales index for food in large retail stores in Spain
#####################################
autoplot(sales, ylab = "Index")
# dev.print(device = pdf, file = "index.pdf")
y <- window(sales, end = c(2018, 12))
yl <- log(y)
u  <- calendar.effects(sales, 
                      trading.day = FALSE, 
                      leap.year = TRUE, 
                      easter = 0)
# Models estimation
m1 <- UC(yl)
m1 <- UC(yl, stepwise = TRUE)
m1 <- UC(yl, stepwise = TRUE, tTest = TRUE)
m2 <- UC(yl, u = u)
m2 <- UC(yl, u = u, stepwise = TRUE)
m2 <- UC(yl, u = u, stepwise = TRUE, tTest = TRUE)
m3 <- UC(yl, u = u, outlier = 4)
m3 <- UC(yl, u = u, outlier = 4, stepwise = TRUE)
m3 <- UC(yl, u = u, outlier = 4, stepwise = TRUE, tTest = TRUE)
# With verbose = TRUE
m1 <- UC(yl, verbose = TRUE)
# Model m3
m3
# Initial parameters
p0 <- getp0(yl, model = "llt/equal/arma(0,0)")
p0[1] <- 0
m <- UC(yl, model = "llt/equal/arma(0,0)", p0 = p0)
# Plotting components
m <- UC(yl, u = window(u, end = c(2018, 12)), outlier = 4)
t <- time(m$comp)
m$comp[is.na(m$comp)] <- 0
p1 <- ggplot() + 
    geom_line(data = yl, aes(x = t, y = yl)) +
    geom_line(data = m$comp[, 1] + m$comp[, 6] + m$comp[, 7], aes(x = t, y = rowSums(m$comp[, c(1, 6, 7)])), color = "red") +
    xlab("") + ylab("Trend")
p2 <- ggplot() + 
    geom_line(data = m$comp[, 3], aes(x = t, y = m$comp[, 3])) +
    xlab("") + ylab("Seasonal")
p3 <- ggplot(data = m$comp[, 4], aes(x = t, y = m$comp[, 4])) + 
    geom_bar(stat = "identity", width = 0.01, fill = "#000000") + 
    xlab("Years") + ylab("Irregular")
grid.arrange(grobs = list(p1, p2, p3),
             widths = 1,
             layout_matrix = rbind(1, 1, 2, 3))
# dev.print(device = pdf, file = "salesComp.pdf")
# Plotting forecasts
m <- UC(yl, u = u, outlier = 4, h = 12)
py <- subset(log(sales), start = 277, end = 300)
forec <- predict(m)
ggplot() + 
  geom_line(data = py, aes(x = time(py), y = py), linetype = "solid") +
  geom_line(data = forec[, 1], aes(x = time(forec), y = forec[, 1]), linetype = "dashed") +
  geom_line(data = forec[, 2], aes(x = time(forec), y = forec[, 2]), linetype = "dotted", color = "red") +
  geom_line(data = forec[, 3], aes(x = time(forec), y = forec[, 3]), linetype = "dotted", color = "red") + 
  geom_vline(xintercept = 2019) +
  xlab("Years") + ylab("Sales (log scale)")
# dev.print(device = pdf, file = "salesFore.pdf")
# Rolling forecasting exercise from observation 228
errors <- rolling(sales, 228)
print(errors$sMAPE)
print(errors$MASE)
# Easter and trading day effect are not statistically significatn strictly speaking
u  <- calendar.effects(sales, 
                       trading.day = TRUE, 
                       leap.year = TRUE, 
                       easter = 06)
m <- UC(yl, u = u, verbose = TRUE)

#####################################
# Example 3: OECD Gross Domestic Product
#            with outliers
#####################################
y <- OECDgdp
# Plot of raw time series
t <- seq(1962, 2019.99, 1 / 4)
ggplot() + 
    geom_line(data = y, aes(x = t, y = y)) + 
    xlab("Years") + ylab("OECD GDP")
# dev.print(device = pdf, file = "gdp.pdf")
yl <- log(y)
# Models
m1 <- UC(yl)
m2 <- UC(yl, model = "?/?/?/?")
m3 <- UC(yl, model = "?/?/arma(2,0)")
m4 <- UC(yl, model = "?/?/arma(2,1)")
# Table 8 in text
firstRow <-     m1$criteria
secondRow <-    m2$criteria
secondRowBis <- m3$criteria
thirdRow <-     m4$criteria
firstRow[4] <-  NA
secondRow[4] <- coef(m2)[5]
ar1 <- -1.5894; ar2 <- 0.6612; rho <- sqrt(ar2); w <- acos(-ar1/2/rho); 
secondRowBis[4] <- 2 * pi / w
ar1 <- -1.5921; ar2 <- 0.6599; rho <- sqrt(ar2); w <- acos(-ar1/2/rho);
thirdRow[4] <- 2 * pi / w
colnames(firstRow) <- c("LLIK", "AIC", "BIC", "Periods")
print(rbind(firstRow, secondRow, secondRowBis, thirdRow))
# Model m2
m2
# Plotting components
comp <- ts(m2$comp[1:232, ], start = start(y), frequency = 4)
p1 <- ggplot() + 
    geom_line(data = y, aes(x = time(comp[, 1]), y = y)) +
    geom_line(data = comp[, 1], aes(x = time(comp[, 1]), y = exp(comp[,1])), color = "red", linetype = "dashed") +
    xlab("") + ylab("Trend")
p2 <- ggplot() + 
    geom_line(data = comp[, 3], aes(x = time(comp[, 3]), y = exp(comp[, 3]))) +
    xlab("Years") + ylab("Cycle")
grid.arrange(p1, p2, nrow = 2)
# dev.print(device = pdf, file = "gdpComp.pdf")
mHP1 <- UChp(yl)
cycles <- ts(exp(cbind(mHP1, m2$comp[1:232, 3], m3$comp[1:232, 3], m4$comp[1:232, 3])), start = start(y), frequency = 4)
ggplot() + 
    geom_line(data = cycles[, 1], aes(x = time(cycles[, 1]), y = cycles[, 1])) +
    geom_line(data = comp[, 3], aes(x = time(comp[, 1]), y = cycles[, 2]), color = "red", linetype = "dashed") +
    geom_line(data = comp[, 3], aes(x = time(comp[, 1]), y = cycles[, 3]), color = "blue", linetype = "dotted") +
    geom_line(data = comp[, 3], aes(x = time(comp[, 1]), y = cycles[, 4]), color = "green", linetype = "twodash") +
    xlab("") + ylab("Trend + Steps")
# Variances and covariances of all cycle estimations
varCycles <- var(cycles)
corCycles <- cor(cycles)
namesVar <- c("HP", "Trigonometric", "AR(2)", "ARMA(2,1)")
colnames(varCycles) <- namesVar
rownames(varCycles) <- namesVar
colnames(corCycles) <- namesVar
rownames(corCycles) <- namesVar
diag(varCycles)
print(corCycles)
# dev.print(device = pdf, file = "gdpCycles.pdf")
#####################################
# Figure 8
#####################################
y1 = UKDriverDeaths
y2 <- ts(c(200.1, 85.2, 180.6, 135.6, 123.3, 159, 127.3, 177, 173.4, 145.7, 135.7, 171.6, 
           144.5, 146.8, 145.2, 109.8, 123.8, 247.8, 83.2, 128.9, 147, 162.8, 145.9, 225.6, 205.8, 148.7, 158.1, 
           156.9, 46.8, 50.3, 59.7, 153.9, 142.3, 124.6, 150.8, 104.7, 130.7, 139.9, 132, 73.6, 78.4, 153.4, 107.7, 
           121.1, 143, 250.5, 249.1, 214.4, 183.9, 86.3, 241.4, 94, 154.5, 87.8, 78.9, 113.6, 118.9, 143, 69.7, 83.4, 
           101.5, 205.1, 137.3, 244.6, 190.5, 151.2, 53, 132.8, 207.7, 131.9, 65.6, 184.7, 249.6, 159.5, 151.3, 184.7, 
           113.7, 157.1, 119.5, 99.5, 123, 110.7, 113.3, 87.9, 93.7, 188.8, 166.1, 82, 131.3, 158.6, 191.1, 144.7, 91.6, 
           78, 104.2, 109, 175, 172.4, 172.6, 138.4, 188.1, 111.4, 74.7, 137.8, 106.8, 103.2, 115.2, 80.6, 122.5, 
           50.4, 149.3, 101.1, 173.7, 125.8, 210.2, 242.8, 163, 128.8, 183.9, 138.5, 180.5, 119.2, 209.3, 129.9, 
           233.1, 251.2, 177.8, 141.7, 194.1, 175.2, 99.6), start = 1849, frequency = 1)
y3 <- ts(c(10, 15, 10, 10, 12, 10, 7, 17, 10, 14, 8, 17, 14, 18, 3, 9, 
           11, 10, 6, 12, 14, 10, 25, 29, 33, 33, 12, 19, 16, 19, 19, 12, 
           34, 15, 36, 29, 26, 21, 17, 19, 13, 20, 24, 12, 6, 14, 6, 12, 9, 
           11, 17, 12, 8, 14, 14, 12, 5, 8, 10, 3, 16, 8, 8, 7, 12, 6, 10, 
           8, 10, 5, 7), frequency = 1)
y4 <- ts(c(116.8, 120.1, 123.2, 130.2, 131.4, 125.6, 124.5,134.3, 135.2, 151.8, 
           146.4, 139, 127.8, 147, 165.9, 165.5, 179.4, 190, 189.8,190.9,203.6, 
           183.5, 169.3, 144.2,141.5, 154.3, 169.5, 193, 203.2, 192.9, 209.4, 227.2, 
           263.7, 297.8, 337.1,361.3, 355.2, 312.6, 309.9, 323.7, 324.1, 355.3, 
           383.4, 395.1, 412.8, 406,438, 446.1, 452.5, 447.3, 475.9, 487.7, 497.2,
           529.8, 551, 581.1, 617.8,658.1, 675.2, 706.6, 724.7), start = 1910, frequency = 1)
y5 <- ts((colSums(matrix(AirPassengers, 3, 144 / 3))),
         start = 1949, frequency = 4)
y6 <- ts(c(303, 169, 152, 257, 247, 189, 146, 220, 248, 195, 141, 235, 278,
           167, 150, 261, 244, 174, 104, 228, 243, 170, 113, 219, 237, 138, 114,
           208, 190, 157, 93, 182, 183, 106, 86, 144, 226, 128, 62, 130, 169, 94,
           91, 188, 148, 114, 62, 139, 104, 99, 76, 122, 107, 76, 51, 111, 95, 71,
           63, 107, 65, 62, 39, 80, 86, 57, 37, 79, 83, 59, 42, 89, 84, 59, 43, 73,
           90, 56, 40, 75, 80, 45, 32, 73, 72, 46, 38, 78, 77, 49, 41, 77, 78, 49,
           34, 72, 68, 51, 23, 42, 72, 49, 39, 64, 63, 43, 45, 56), start = 1960, frequency = 4)
p1 <- ggplot() + 
    geom_line(data = y1, aes(x = time(y1), y = y1))+
    xlab("") + ylab("UK Deaths")
p2 <- ggplot() + 
    geom_line(data = y2, aes(x = time(y2), y = y2))+
    xlab("") + ylab("Rainfall in Brazil")
p3 <- ggplot() + 
    geom_line(data = y3, aes(x = time(y3), y = y3))+
    xlab("") + ylab("Purse snatching")
p4 <- ggplot() + 
    geom_line(data = y4, aes(x = time(y4), y = y4))+
    xlab("") + ylab("US GNP")
p5 <- ggplot() + 
    geom_line(data = y5, aes(x = time(y5), y = y5))+
    xlab("") + ylab("Air passengers")
p6 <- ggplot() + 
    geom_line(data = y6, aes(x = time(y6), y = y6))+
    xlab("") + ylab("Coal consumption")
grid.arrange(grobs = list(p1, p2, p3, p4, p5, p6),
             layout_matrix = rbind(c(1, 2), c(3, 4), c(5, 6)))
# dev.print(device=pdf, file="figure8.pdf")
#####################################
# Examples from other papers / books
#####################################
# UK drivers killed and seriously injured
# Harvey (1989) pp. 82
y <- log(UKDriverDeaths)
# Model in Harvey (1989), pp. 83
m1 <- UC(subset(y, end = 168), model = "llt/equal/arma(0,0)")
# Automatic identification
m2 <- UC(subset(y, end = 168))
# Model with outliers and full sample
m3 <- UC(y, outlier = 4, model = "rw/equal/arma(0,0)")
# Optimal model with outliers
m4 <- UC(y, outlier = 4)
# Table 9
table9 <- matrix(NA, 8, 4)
colnames(table9) <- c("STAMP & UCOMP", "UComp optimal", "UComp", "UComp optimal")
rownames(table9) <- c("Level", "Slope", "Seasonal", "Irregular", "LS Feb1983", "LLIK", "AIC", "BIC")
table9[c(1,2,3,4,6,7,8), 1] <-  c(coef(m1), m1$criteria[1:3])
table9[c(1,3,4,6,7,8), 2] <-  c(coef(m2), m2$criteria[1:3])
table9[c(1,3,4,5,6,7,8), 3] <-  c(coef(m3), m3$criteria[1:3])
table9[c(1,3,4,5,6,7,8), 4] <-  c(coef(m4), m4$criteria[1:3])
print(table9)
## Rainfall in Fortaleza
#  Harvey (1989), pp. 86
y <- ts(c(200.1, 85.2, 180.6, 135.6, 123.3, 159, 127.3, 177, 173.4, 145.7, 135.7, 171.6, 
          144.5, 146.8, 145.2, 109.8, 123.8, 247.8, 83.2, 128.9, 147, 162.8, 145.9, 225.6, 205.8, 148.7, 158.1, 
          156.9, 46.8, 50.3, 59.7, 153.9, 142.3, 124.6, 150.8, 104.7, 130.7, 139.9, 132, 73.6, 78.4, 153.4, 107.7, 
          121.1, 143, 250.5, 249.1, 214.4, 183.9, 86.3, 241.4, 94, 154.5, 87.8, 78.9, 113.6, 118.9, 143, 69.7, 83.4, 
          101.5, 205.1, 137.3, 244.6, 190.5, 151.2, 53, 132.8, 207.7, 131.9, 65.6, 184.7, 249.6, 159.5, 151.3, 184.7, 
          113.7, 157.1, 119.5, 99.5, 123, 110.7, 113.3, 87.9, 93.7, 188.8, 166.1, 82, 131.3, 158.6, 191.1, 144.7, 91.6, 
          78, 104.2, 109, 175, 172.4, 172.6, 138.4, 188.1, 111.4, 74.7, 137.8, 106.8, 103.2, 115.2, 80.6, 122.5, 
          50.4, 149.3, 101.1, 173.7, 125.8, 210.2, 242.8, 163, 128.8, 183.9, 138.5, 180.5, 119.2, 209.3, 129.9, 
          233.1, 251.2, 177.8, 141.7, 194.1, 175.2, 99.6), start = 1849, frequency = 1)
m1 <- UC(y, model = "?/-8/?/?")
m2 <- UC(y, model = "?/?/?/?")
# Left half of Table 10
table10l <- matrix(NA, 9, 2)
colnames(table10l) <- c("STAMP & UCOMP", "UComp optimal")
rownames(table10l) <- c("Rho", "Period", "Var", "Irregular", "AR(1)", "Const.", "LLIK", "AIC", "BIC")
table10l[6, 1] <- m1$comp[1]
table10l[c(1,2,3,4,7,8,9), 1] <-  c(coef(m1), m1$criteria[1:3])
table10l[c(4,5,6,7,8,9), 2] <-  c(coef(m2), m2$criteria[1:3])
print(table10l)
## Purse snatchings in Hyde Park, Chicago
#  Harvey (1989), pp. 89
y <- ts(c(10, 15, 10, 10, 12, 10, 7, 17, 10, 14, 8, 17, 14, 18, 3, 9, 
          11, 10, 6, 12, 14, 10, 25, 29, 33, 33, 12, 19, 16, 19, 19, 12, 
          34, 15, 36, 29, 26, 21, 17, 19, 13, 20, 24, 12, 6, 14, 6, 12, 9, 
          11, 17, 12, 8, 14, 14, 12, 5, 8, 10, 3, 16, 8, 8, 7, 12, 6, 10, 
          8, 10, 5, 7), frequency = 1)
m1 <- UC(y, model = "rw/none/arma(0,0)")
m2 <- UC(y)
# Right half of Table 10
table10r <- matrix(NA, 8, 2)
colnames(table10r) <- c("STAMP & UCOMP", "UComp optimal")
rownames(table10r) <- c("Level", "AR(1)", "AR(2)", "Irregular", "Const.", "LLIK", "AIC", "BIC")
table10r[c(1,4,6,7,8), 1] <-  c(coef(m1), m1$criteria[1:3])
table10r[c(4,2,3,5,6,7,8), 2] <-  c(coef(m2), m2$criteria[1:3])
print(table10r)
## US GNP data
# Harvey (1989), pp. 92
y <- ts(c(116.8, 120.1, 123.2, 130.2, 131.4, 125.6, 124.5,134.3, 135.2, 151.8, 
          146.4, 139, 127.8, 147, 165.9, 165.5, 179.4, 190, 189.8,190.9,203.6, 
          183.5, 169.3, 144.2,141.5, 154.3, 169.5, 193, 203.2, 192.9, 209.4, 227.2, 
          263.7, 297.8, 337.1,361.3, 355.2, 312.6, 309.9, 323.7, 324.1, 355.3, 
          383.4, 395.1, 412.8, 406,438, 446.1, 452.5, 447.3, 475.9, 487.7, 497.2,
          529.8, 551, 581.1, 617.8,658.1, 675.2, 706.6, 724.7), start = 1910, frequency = 1)
m1 <- UC(log(y), model = "llt/-8/none/arma(0,0)")
m2 <- UC(log(y), model = "?/?/?/?")
# Left half of Table 11
table11l <- matrix(NA, 10, 2)
colnames(table11l) <- c("STAMP & UCOMP", "UComp optimal")
rownames(table11l) <- c("Damping", "Level", "Slope", "Rho", "Period", "Var", "Irregular", "LLIK", "AIC", "BIC")
table11l[c(2:10), 1] <-  c(coef(m1), m1$criteria[1:3])
table11l[c(1:6,8:10), 2] <-  c(coef(m2), m2$criteria[1:3])
print(table11l)
## Quartely air passengers data
# Harvey (1989), pp. 93
y <- ts(log(colSums(matrix(AirPassengers, 3, 144 / 3))),
        start = 1949, frequency = 4)
m1 <- UC(y, model = "llt/equal/arma(0,0)")
m2 <- UC(y)
# Rigt half of Table 11
table11r <- matrix(NA, 9, 2)
colnames(table11r) <- c("STAMP & UCOMP", "UComp optimal")
rownames(table11r) <- c("Level", "Slope", "Seasonal", "Season(4)", "Season(2)", "Irregular", "LLIK", "AIC", "BIC")
table11r[c(1,2,3,6:9), 1] <-  c(coef(m1), m1$criteria[1:3])
table11r[c(1,2,4,5,7,8,9), 2] <-  c(coef(m2), m2$criteria[1:3])
print(table11r)
## UK coal consumption in the UK
# Harvey (1989), pp. 95
y <- ts(c(303, 169, 152, 257, 247, 189, 146, 220, 248, 195, 141, 235, 278,
          167, 150, 261, 244, 174, 104, 228, 243, 170, 113, 219, 237, 138, 114,
          208, 190, 157, 93, 182, 183, 106, 86, 144, 226, 128, 62, 130, 169, 94,
          91, 188, 148, 114, 62, 139, 104, 99, 76, 122, 107, 76, 51, 111, 95, 71,
          63, 107, 65, 62, 39, 80, 86, 57, 37, 79, 83, 59, 42, 89, 84, 59, 43, 73,
          90, 56, 40, 75, 80, 45, 32, 73, 72, 46, 38, 78, 77, 49, 41, 77, 78, 49,
          34, 72, 68, 51, 23, 42, 72, 49, 39, 64, 63, 43, 45, 56), start = 1960, frequency = 4)
m1 <- UC(log(subset(y, end = 92)), model = "llt/eq/arma(0,0)")
m2 <- UC(log(subset(y, end = 92)))
# Right half of Table 11
table12 <- matrix(NA, 7, 2)
colnames(table12) <- c("STAMP & UCOMP", "UComp optimal")
rownames(table12) <- c("Level", "Slope", "Seasonal", "Irregular", "LLIK", "AIC", "BIC")
table12[, 1] <-  c(coef(m1), m1$criteria[1:3])
table12[c(1,3:7), 2] <-  c(coef(m2), m2$criteria[1:3])
print(table12)