dtu02417a3: ARIMAX Model for Building Data

Question 3.1

frameData <- function(){
    dat <- read.csv("house_data_30min.csv")
    time <- as.numeric(dat[ , 1])
    heating <- as.numeric(dat[ , 2])
    tempExternal <- as.numeric(dat[ , 3])
    iSolar <- as.numeric(dat[ , 4])
    series <- seq(1, length(time))
    dat.f <- data.frame(series, time, heating, tempExternal, iSolar)
    return(dat.f)
}
dat.f <- frameData()
rm(frameData)
n <- length(dat.f$series) - 4

Read the data and plot the given quantities as a function of time. Indicate the test data in the plots. Comment on the evolution of the values over time.

plot(dat.f$series, dat.f$heating, type = "l", col = "blue", lwd = 1, lty = 1, cex = 0.8, bty = "n", 
     main = "Training and Testing Data of Heating Power Inside the Building", 
     xlab = "Series", ylab = "Heat Power Inside (W)")
points(dat.f$series[(n+1):(n+4)], dat.f$heating[(n+1):(n+4)], 
       col = "blue", lty = 1, pch = 16, cex = 0.8)

legend("topleft", inset = .02, legend = c("Training Data", "Testing Data"), col = c("blue", "blue"), 
       pch = c(NaN, 16), lty = c(1, 1), lwd = c(2, 2))

plot(dat.f$series, dat.f$tempExternal, type = "l", col = "blue", lwd = 1, lty = 1, cex = 0.8, bty = "n", 
     main = "Training and Testing Data of External Temperature", xlab = "Series", ylab = "Temperature (Celsius D)")
points(dat.f$series[(n+1):(n+8)], dat.f$tempExternal[(n+1):(n+8)], 
       col = "blue", lty = 1, pch = 16, cex = 0.8)

legend("bottomleft", inset = .02, legend = c("Training Data", "Testing Data"), col = c("blue", "blue"), 
       pch = c(NaN, 16), lty = c(1, 1), lwd = c(2, 2))

plot(dat.f$series, dat.f$iSolar, type = "l", col = "blue", lwd = 1, lty = 1, cex = 0.8, bty = "n", 
     main = "Training and Testing Data of Solar Irradiation", xlab = "Series", ylab = "Solar Irradiation (W / m2)")
points(dat.f$series[(n+1):(n+8)], dat.f$iSolar[(n+1):(n+8)], 
       col = "blue", lty = 1, pch = 16, cex = 0.8)

legend("topleft", inset = .02, legend = c("Training Data", "Testing Data"), col = c("blue", "blue"), 
       pch = c(NaN, 16), lty = c(1, 1), lwd = c(2, 2))

Question 3.2

# Plot the time series / residuals and its ACF, PACF, and Ljung-Box plot
plotTimeSeriesResidual <- function(dat, num_lag = 20, ...){
    # If the input dat is arima model, just take its residuals are the dat.
    if(class(dat) == "Arima")
        dat <- dat$residuals
    par(mfrow = c(4, 1), mgp = c(2, 0.7, 0), mar = c(1, 1, 1, 1), oma = c(2, 2, 4, 2), 
        cex.lab = 0.8, cex.axis = 0.8)
    plot(dat, type = "l", bty = "n")
    title(main = paste("TS, ACF, PACF and p Values for Ljung-Box Statistic"), cex.main = 1.5, line = 1, outer = TRUE)
    # ACF and PACF
    acf(dat, lag.max = num_lag, xaxt = "n", bty = "n")
    axis(1, at = seq(0, 20, by = 1), las = 0, outer = FALSE)
    pacf(dat, lag.max = num_lag, xaxt = "n", bty = "n")
    axis(1, at = seq(1, 20, by = 1), las = 0, outer = FALSE)
    # Ljung-Box Plot
    value_p <- sapply(1: num_lag, function(i) Box.test(dat, i, type = "Ljung-Box")$p.value)
    plot(1L: num_lag, value_p, xlab = "lag", ylab = "p value", ylim = c(0, 1), bty = "n", xaxt = "n")
    axis(1, at = seq(1, 20, by = 1), las = 0, outer = FALSE)
    abline(h = 0.05, lty = 2, col = "blue")
}

plotTimeSeriesResidual(dat.f$heating)

Question 3.3: ARIMA Model

mod_arima_1 <- arima(dat.f$heating, order = c(1, 0, 0))
mod_arima_1

Call:
arima(x = dat.f$heating, order = c(1, 0, 0))

Coefficients:
         ar1  intercept
      0.9002    73.5528
s.e.  0.0195     3.5272

sigma^2 estimated as 62.43:  log likelihood = -1695,  aic = 3395.99

plotTimeSeriesResidual(mod_arima_1)

mod_arima_2 <- arima(dat.f$heating, order = c(1, 0, 0), include.mean = F)
mod_arima_2

Call:
arima(x = dat.f$heating, order = c(1, 0, 0), include.mean = F)

Coefficients:
         ar1
      0.9938
s.e.  0.0042

sigma^2 estimated as 65.52:  log likelihood = -1708.11,  aic = 3420.23

plotTimeSeriesResidual(mod_arima_2)

mod_arima_3 <- arima(dat.f$heating, order = c(1, 0, 2))
mod_arima_3

Call:
arima(x = dat.f$heating, order = c(1, 0, 2))

Coefficients:
         ar1     ma1      ma2  intercept
      0.9037  0.0742  -0.0932    73.5135
s.e.  0.0247  0.0544   0.0503     3.5590

sigma^2 estimated as 61.42:  log likelihood = -1691.06,  aic = 3392.12

plotTimeSeriesResidual(mod_arima_3)