dtu02417a3: ARIMA Model for Building Data

Question 3.3, ARIMA Model of Heating

mod_heating <- arima(dat.f$heating[1: num_obs], order = c(0, 1, 2), seasonal = list(order = c(0, 0, 2), period = 48), 
                     fixed = c(0, NA, NA, NA), include.mean = T)
mod_heating

Call:
arima(x = dat.f$heating[1:num_obs], order = c(0, 1, 2), seasonal = list(order = c(0, 
    0, 2), period = 48), include.mean = T, fixed = c(0, NA, NA, NA))

Coefficients:
      ma1      ma2    sma1    sma2
        0  -0.1261  0.0696  0.1217
s.e.    0   0.0490  0.0479  0.0474

sigma^2 estimated as 56.66:  log likelihood = -1654.23,  aic = 3316.47

plotTimeSeriesResidual(mod_heating, num_lag = 100, "ARIMA Model 3.1.4 of Heating")

setEPS()
postscript("7.eps", width = 11, height = 12)
plotTimeSeriesResidual(mod_heating, num_lag = 100, "ARIMA Model 3.1.4 of Heating")
dev.off()

pred_heating <- predict(mod_heating, n.ahead = 4)
intervalPred <- calIntervalPred(482, 13, pred_heating$pred, pred_heating$se)
mat_intervalPred <- matrix(nrow = 2, ncol = 4)
mat_intervalPred[1,] <- intervalPred$boundUp
mat_intervalPred[2,] <- intervalPred$boundLow
plotPredHeating(pred_heating$pred, mat_intervalPred, str = "ARIMA(0,1,2)(0,0,1)48 Model")

rm(mat_intervalPred)

setEPS()
postscript("12.eps", width = 11, height = 6)
plotPredHeating(pred_heating$pred, mat_intervalPred, str = "ARIMA(0,1,2)(0,0,1)48 Model")
dev.off()
rm(mat_intervalPred)

Question 3.4, ARIMAX Model of Heating

ts_tempExternal <- ts(dat.f$tempExternal[1: num_obs], frequency = 48)
ts_iSolar <- ts(dat.f$iSolar[1: num_obs], frequency = 48)
ts_heating <- ts(dat.f$heating[1: num_obs], frequency = 48)
auto.arima(ts_tempExternal, d = 1, seasonal = TRUE)
auto.arima(ts_iSolar, d = 0, seasonal = TRUE)
auto.arima(ts_heating, xreg = ts_tempExternal, d = 1, seasonal = TRUE)
auto.arima(ts_heating, xreg = ts_iSolar, d = 0, seasonal = TRUE)

Pre-Whitening of Heating using ARIMA Model of External Temperature

mod_tempExternal <- arima(dat.f$tempExternal[1: num_obs], order = c(0, 1, 1), seasonal = list(order = c(0, 0, 1), 
    period = 48), include.mean = TRUE)
x_1 <- dat.f$tempExternal[1: num_obs] - fitted(Arima(dat.f$tempExternal[1: num_obs], model = mod_tempExternal))
y_1 <- dat.f$heating[1: num_obs] - fitted(Arima(dat.f$heating[1: num_obs], model = mod_tempExternal))
sum(x_1 - mod_tempExternal$residuals) < 10^6

[1] TRUE

par(mar=c(3.5,3.5,3,1), mgp=c(2,0.7,0))
plotCCFunc(x_1, y_1, "Filtered External Temp", "Filtered Heating")

setEPS()
postscript("8.eps", width = 11, height = 12)
plotTimeSeriesResidual(mod_tempExternal, num_lag = 100, "ARIMA Model 1.1.1 of External Temperature")
dev.off()

setEPS()
postscript("10.eps", width = 11, height = 6)
plotCCFunc(x_1, y_1, "Filtered External Temp", "Filtered Heating")
dev.off()

Pre-Whitening of Heating using ARIMA Model of Solar Irradiation

mod_iSolar <- arima(dat.f$iSolar[1: num_obs], order = c(4, 0, 0), seasonal = list(order = c(0, 0, 2), period = 48), 
                    include.mean = T)
x_2 <- dat.f$iSolar[1: num_obs] - fitted(Arima(dat.f$iSolar[1: num_obs], model = mod_iSolar))
y_2 <- residuals(Arima(dat.f$heating[1: num_obs], model = mod_iSolar))
sum(x_2 - mod_tempExternal$residuals) < 10^6

[1] TRUE

par(mar=c(3.5,3.5,3,1), mgp=c(2,0.7,0))
plotCCFunc(x_2, y_2, "Filtered Solar Irradiation", "Filtered Heating")

setEPS()
postscript("11.eps", width = 11, height = 6)
plotCCFunc(x_2, y_2, "Filtered Solar Irradiation", "Filtered Heating")
dev.off()

Prediction of Heating using ARIMAX Model

ts_heating <- ts(dat.f$heating[1: num_obs], frequency = 48)
ts_iSolar <- ts(dat.f$iSolar[1: num_obs], frequency = 48)
mat_obs <- cbind(y = ts_heating, x_1 = ts_iSolar, x_2 = lag(ts_iSolar, k = - 1))
mod_heating_x <- lm(y ~ x_1 + x_2, data = mat_obs, na.action = na.omit)

mod_residualHeating <- auto.arima(residualHeating[1: 481])
pred_residualHeating <- predict(mod_residualHeating, n.ahead = 4)$pred

mat_x <- cbind(rep(1, 4), dat.f$iSolar[(num_obs+1): (num_obs+4)], dat.f$iSolar[(num_obs): (num_obs+3)])
pred_heating_x <- mat_x %*% mod_heating_x$coefficients
pred_heating_x_residual <- pred_heating_x + pred_residualHeating
# plotPredHeating(pred_heating_x, "ARIMAX Model with Known Solar Irradiation as Input")
plotPredHeating(pred_heating_x_residual, "ARIMAX Model with Known Solar Irradiation as Input")

setEPS()
postscript("13.eps", width = 11, height = 6)
plotPredHeating(pred_heating_x_residual, "ARIMAX Model with Known Solar Irradiation as Input")
dev.off()

mod_heating_3

Call:
arima(x = dat.f$heating[2:num_obs], order = c(3, 0, 4), seasonal = list(order = c(0, 
    0, 1), period = 48), xreg = cbind(tempExternal = dat.f$tempExternal[2:num_obs], 
    iSolar = dat.f$iSolar[2:num_obs], iSolar_1 = dat.f$iSolar[1:num_obs - 1]))

Coefficients:
          ar1     ar2     ar3     ma1     ma2      ma3     ma4    sma1  intercept  tempExternal  iSolar  iSolar_1
      -0.3056  0.2911  0.8903  0.6817  0.3588  -0.2915  0.2071  0.1442    87.3190       -1.6666  -0.067   -0.0368
s.e.   0.0989  0.0668  0.0776  0.1064  0.0780   0.0720  0.0500  0.0410     3.3126        0.4677   0.002    0.0019

sigma^2 estimated as 10.73:  log likelihood = -1254.8,  aic = 2535.6

pred_heating_3 <- predict(mod_heating_3, newxreg = cbind(tempExternal = 
    dat.f$tempExternal[(num_obs + 1): (num_obs + 4)], iSolar = dat.f$iSolar[(num_obs + 1): (num_obs + 4)], 
    iSolar_1 = dat.f$iSolar[num_obs: (num_obs + 3)]), n.ahead = 4)$pred
plotPredHeating(pred_heating_3, "ARIMAX Model with Known Solar Irradiation as Input")

ARIMAX Model using External Temperature and Solar Irradiation as Input

mod_heating_4

Call:
arima(x = dat.f$heating[2:num_obs], order = c(3, 0, 4), seasonal = list(order = c(0, 
    0, 2), period = 48), xreg = cbind(tempExternal = dat.f$tempExternal[2:num_obs], 
    iSolar = dat.f$iSolar[2:num_obs], iSolar_1 = dat.f$iSolar[1:num_obs - 1]))

Coefficients:
          ar1     ar2     ar3     ma1     ma2      ma3     ma4    sma1    sma2  intercept  tempExternal   iSolar
      -0.3030  0.2717  0.9089  0.6855  0.3925  -0.2983  0.1985  0.1568  0.2190    86.8396       -1.5537  -0.0677
s.e.   0.1117  0.1127  0.0751  0.1198  0.0998   0.0728  0.0500  0.0477  0.0511     3.5250        0.4309   0.0019
      iSolar_1
       -0.0373
s.e.    0.0018

sigma^2 estimated as 10.23:  log likelihood = -1245.55,  aic = 2519.1

setEPS()
postscript("15.eps", width = 11, height = 12)
plotTimeSeriesResidual(mod_heating_4$residuals, num_lag = 100)
dev.off()

plotPredHeating(pred_heating_4$pred, mat_intervalPred, str = "ARIMAX Model with Known Ex-Temp and Solar-Irra as Input")

There were 25 warnings (use warnings() to see them)

setEPS()
postscript("14.eps", width = 11, height = 6)
plotPredHeating(pred_heating_4$pred, mat_intervalPred, str = "ARIMAX Model with Known Ex-Temp and Solar-Irra as Input")
dev.off()

mod_heating_5

Call:
arima(x = dat.f$heating[3:num_obs], order = c(3, 0, 4), seasonal = list(order = c(0, 
    0, 1), period = 48), xreg = cbind(tempExternal = dat.f$tempExternal[3:num_obs], 
    iSolar = dat.f$iSolar[3:num_obs], iSolar_1 = dat.f$iSolar[2:(num_obs - 1)], 
    iSolar_2 = dat.f$iSolar[1:(num_obs - 2)]))

Coefficients:
          ar1     ar2     ar3     ma1     ma2      ma3     ma4    sma1  intercept  tempExternal   iSolar  iSolar_1
      -0.3681  0.2881  0.9419  0.7496  0.4042  -0.3146  0.2103  0.1534    87.5979       -1.7410  -0.0664   -0.0375
s.e.   0.0247  0.0322  0.0237  0.0702  0.0785   0.0623  0.0529  0.0405     3.2543        0.4851   0.0021    0.0020
      iSolar_2
        0.0021
s.e.    0.0019

sigma^2 estimated as 10.65:  log likelihood = -1251.88,  aic = 2531.76

pred_heating_6 <- predict(mod_heating_6, newxreg = cbind(tempExternal = 
    dat.f$tempExternal[(num_obs + 1): (num_obs + 4)], iSolar = dat.f$iSolar[(num_obs + 1): (num_obs + 4)], 
    iSolar_1 = dat.f$iSolar[num_obs: (num_obs + 3)]), n.ahead = 4)$pred
plotPredHeating(pred_heating_6, "ARIMAX Model with Known Solar Irradiation as Input")