R Notebook

library(datasets)
library(fUnitRoots)
library(forecast)
library(astsa)
?UKgas
ukgas_ts <- UKgas
class(ukgas_ts)

[1] "ts"

No need to convert to time series– it already is in time series format “ts”.

Preliminary Discussion and Analysis

First, we will plot the uk_gas time series data and observe trends

par(mfrow = c(1,1), cex.axis = 0.9)

plot(ukgas_ts,
     main = "Gas Consumption in the UK (Quarterly, 1960-1986)",
     xlab = 'Year', ylab = 'Millions of Therms (Natural Gas)',
     cex.axis = 0.5)

axis(1, at = time(ukgas_ts)[cycle(ukgas_ts) == 1], labels = FALSE) #has marks for year

Inititial Plot an Brief Discussion: Quarterly time series data From Q1 1960 to Q4 1986 Natural gas consumption is measured in millions of therms

Variance is increasing as mean increases. from 1960 to 1967, variance appears roughly the same. Then it seems to increase exponentially once mean starts to increase. Seaonality is a key factor here. Perhaps more people using natural gas to heat their homes, starting in the 70s.

UK National project to convert from town gas (made from coal) to natural gas between 1966 and 1977. In 1970 a huge supply of Natural gas discovered in North Sea.

This change in variance, along with trend and seasonality, may make this fairly complex to model.

There is an upward trend. Seemingly flat from 1960-1066, but then increaing in future years.

A variance transformation is absolutely necessary here

First, we will try to stabilize variance with a log transformation

# Take log transformation in attempts to stabilize variance
log_ukgas_ts <- log(ukgas_ts)
plot(log_ukgas_ts, main = "Log ukgas Prices", ylab = "Log Price") 

axis(1, at = time(ukgas_ts)[cycle(ukgas_ts) == 1], labels = FALSE) #has marks for year

Variance doesnt appear to be fully stabilized with a log transformation.

So we will do a box cox transformation instead.

Estimate value for lambda. This function finds the best lambda.

lambda <- BoxCox.lambda(ukgas_ts)
lambda

[1] -0.4457023

Apply a box cox transformation (using best lambda) to the data

bc_ukgas_ts <- BoxCox(ukgas_ts, lambda = lambda)
plot(bc_ukgas_ts,
     main = "Box-Cox Transformed UKgas Time Series (lambda = -0.445)", 
     ylab = "Transformed gas consumption",
     xlab = "Year")

# Add tick marks for every year
axis(1, at = time(bc_ukgas_ts)[cycle(bc_ukgas_ts) == 1], labels = FALSE)

This looks much better than the log transformation. Variance appears to be much more stabilized.

However, we still have seasonality and trend that must be accounted for to achieve stationarity.

We can test for stationarity about the mean with the ADF Test.

Its obviously not stationary right now; there’s a clear trend up. But let’s do it anyway.

# Get the number of observations in the TS
(length(ukgas_ts)-1)^(1/3)

[1] 4.747459

So we will round down and go with lag=4 for adfTest()

library(fUnitRoots)
adfTest(bc_ukgas_ts, lags = 4, type = 'c')

Title:
 Augmented Dickey-Fuller Test

Test Results:
  PARAMETER:
    Lag Order: 4
  STATISTIC:
    Dickey-Fuller: -0.4562
  P VALUE:
    0.8866 

Description:
 Wed Jul 23 17:46:09 2025 by user: Sferg

With p-value > 0.05,

Not stationary. Must difference again. Trend is still there.

There is both seasonality and trend. To account for quarterly seasonality, a lag 4 difference will be applied. Although the primary purporse of taking a lag 4 difference is to remove seasonality, the lag 4 difference might be enough to also take care of the trend.

First, I’ll take the lag 4 difference. Then, I’ll use the ADF test to determine if we must difference again to remove trend.

In efforts to be parsimonious, it is ideal to only take one difference if possible.

# Seasonal difference only
diff4_bc_ukgas_ts <- diff(bc_ukgas_ts, lag=4)
plot(diff4_bc_ukgas_ts,
     main = "Seasonally Differenced, Box-Cox Transformed \n UKgas Time Series (lambda = -0.445)", 
     ylab = "Transformed gas consumption",
     xlab = "Year")

# Line at y=0
abline(h = 0, col = "red", lwd = .5)

# Add tick marks for every year
axis(1, at = time(bc_ukgas_ts)[cycle(bc_ukgas_ts) == 1], labels = FALSE)

Note: Spike in volatility post-1970. Outlier period 1970-1972. In 1970 BP discovered a huge supply of natural gas in the North Sea. This volatility could be a consequence of that, combined with developing infrastructure. ” While the spike around 1970–1972 aligns more with North Sea gas development, any lingering post-1973 anomalies could reflect the broader 1973–1974 energy crisis, including demand shocks, hoarding behavior, or changes in industrial output. ”

All in all, weird stuff was happening in the UK natural gas market around 1970.

adfTest(diff4_bc_ukgas_ts, lags=4, type = 'nc')

Title:
 Augmented Dickey-Fuller Test

Test Results:
  PARAMETER:
    Lag Order: 4
  STATISTIC:
    Dickey-Fuller: -2.0099
  P VALUE:
    0.04484 

Description:
 Wed Jul 23 17:46:09 2025 by user: Sferg

Using a conservative ADF specification with no constant, the p-value was < 0.05, indicating that the Box-Cox transformed and lag-4 differenced UKgas series is stationary.”

The ADF test result narrowly rejects the null hypothesis of a unit root at the 0.05 significance level. This suggests that the Box-Cox transformation (\(\lambda\) = –0.4457) and seasonal lag-4 differencing were sufficient to stabilize the mean and achieve stationarity — at least according to the test.

That said, the result is borderline, and some residual trend may still be present. We’ll note this as a possible limitation, but proceed with modeling under the assumption of weak stationarity based on this result.

Note, we could apply a regular difference to the data that has already been lag4 difference and bc transformed. If we did that, the p-value would shrink, and there would be no doubt that this would make our data stationary about the mean.

However, in efforts to avoid conducting any unnecessary transformations, we lean on the results of the ADF test. The p-value is 0.0448. Although that is close to 0.05, we will conclude that we conducted sufficient transformations to conclude stationarity.

In order to speculate on the ARIMA structure of the time series data, we’ll look at the ACF and PACF plots of our transformed data.

#Look at ACF and PACF of the lag-4 differenced data
par(mfrow = c(1,2))
acf(diff4_bc_ukgas_ts, lag.max = 20, xlab = "LAG ÷ 4") # 4 years
pacf(diff4_bc_ukgas_ts, lag.max = 20, xlab = "LAG ÷ 4")

Each unit on the x-axis represents 1 year, since the data is quarterly and the seasonal period is 4.

In athe ACF plot, Lags 1, 4, and 5 show significant autocorrelations. While seasonal differencings (of lag=4) removed much of the seasonality, some residual seasonal structure reamains.

Lag 1 spike suggests short-term autocorellation : Perhaps MA(1) or AR(1) With a spike again at Lag 4, this suggests some seasonal autocorrelation. So maybe an SAR(1) or SMA(1) term

The PACF has signicant spikes at quarters 1, 4, 5, 9

Therefore, there is some residual seasonal structure even after the seasonal (lag4) differencing. However, the remaining lags fall mostly within the 95% CIs. This suggests weak/no autocorrelation after 1 year (4 lags).

So what we did: 1. Box Cox transformation with lamda = -0.445 2. Seasonal differencing, lag = 4

Based on my model transformations and the resulting ACF and PACF plots, we can speculate the structure of the models that auto.arima() will generate for us.

We applied a seasonal difference of lag 4. Therefore, D=1 and the lag term [4] We did not apply a regular difference. Therefore, d=1.

Make some predictions of what our model might look like

We now use auto.arima on the Box-Cox transformed data. We decided to feed auto.arima() the un-differenced data. The auto-arima function should handle seasonal differencing on its own (with a D=1 term).

We are allowing drift and we display the trace of auto.arima() to evaluate other potentially suitable models with similar AICc values.

After reviewing the documenation for auto.arima(), we decided to add the command “stepwise=FALSE”: giving us a much more comprehensive model search.

We will evaluate 4 models. 4 Arima models before choosing 1 to accept and use for forecasting.

# Apply the 
auto.arima(bc_ukgas_ts, trace = TRUE, approximation = FALSE, allowdrift = TRUE, stepwise=FALSE)

 ARIMA(0,0,0)(0,1,0)[4]                    : -637.6961
 ARIMA(0,0,0)(0,1,0)[4] with drift         : -663.1634
 ARIMA(0,0,0)(0,1,1)[4]                    : -635.6722
 ARIMA(0,0,0)(0,1,1)[4] with drift         : -668.4061
 ARIMA(0,0,0)(0,1,2)[4]                    : -639.1875
 ARIMA(0,0,0)(0,1,2)[4] with drift         : -666.6154
 ARIMA(0,0,0)(1,1,0)[4]                    : -635.7008
 ARIMA(0,0,0)(1,1,0)[4] with drift         : -668.7646
 ARIMA(0,0,0)(1,1,1)[4]                    : Inf
 ARIMA(0,0,0)(1,1,1)[4] with drift         : -666.6188
 ARIMA(0,0,0)(1,1,2)[4]                    : Inf
 ARIMA(0,0,0)(1,1,2)[4] with drift         : -664.4158
 ARIMA(0,0,0)(2,1,0)[4]                    : -641.0757
 ARIMA(0,0,0)(2,1,0)[4] with drift         : -666.6203
 ARIMA(0,0,0)(2,1,1)[4]                    : Inf
 ARIMA(0,0,0)(2,1,1)[4] with drift         : Inf
 ARIMA(0,0,0)(2,1,2)[4]                    : Inf
 ARIMA(0,0,0)(2,1,2)[4] with drift         : Inf
 ARIMA(0,0,1)(0,1,0)[4]                    : -636.0145
 ARIMA(0,0,1)(0,1,0)[4] with drift         : -665.4598
 ARIMA(0,0,1)(0,1,1)[4]                    : -633.9025
 ARIMA(0,0,1)(0,1,1)[4] with drift         : -667.8465
 ARIMA(0,0,1)(0,1,2)[4]                    : -637.0249
 ARIMA(0,0,1)(0,1,2)[4] with drift         : -666.3497
 ARIMA(0,0,1)(1,1,0)[4]                    : -633.9094
 ARIMA(0,0,1)(1,1,0)[4] with drift         : -668.4184
 ARIMA(0,0,1)(1,1,1)[4]                    : -632.3564
 ARIMA(0,0,1)(1,1,1)[4] with drift         : -666.2778
 ARIMA(0,0,1)(1,1,2)[4]                    : -650.412
 ARIMA(0,0,1)(1,1,2)[4] with drift         : -664.0973
 ARIMA(0,0,1)(2,1,0)[4]                    : -639.7273
 ARIMA(0,0,1)(2,1,0)[4] with drift         : -666.2982
 ARIMA(0,0,1)(2,1,1)[4]                    : Inf
 ARIMA(0,0,1)(2,1,1)[4] with drift         : -664.0747
 ARIMA(0,0,1)(2,1,2)[4]                    : Inf
 ARIMA(0,0,1)(2,1,2)[4] with drift         : Inf
 ARIMA(0,0,2)(0,1,0)[4]                    : -641.0514
 ARIMA(0,0,2)(0,1,0)[4] with drift         : -663.6582
 ARIMA(0,0,2)(0,1,1)[4]                    : -638.8933
 ARIMA(0,0,2)(0,1,1)[4] with drift         : -665.8615
 ARIMA(0,0,2)(0,1,2)[4]                    : -639.5629
 ARIMA(0,0,2)(0,1,2)[4] with drift         : -664.1997
 ARIMA(0,0,2)(1,1,0)[4]                    : -638.8963
 ARIMA(0,0,2)(1,1,0)[4] with drift         : -666.3627
 ARIMA(0,0,2)(1,1,1)[4]                    : -636.9265
 ARIMA(0,0,2)(1,1,1)[4] with drift         : -664.1691
 ARIMA(0,0,2)(1,1,2)[4]                    : Inf
 ARIMA(0,0,2)(1,1,2)[4] with drift         : -661.9109
 ARIMA(0,0,2)(2,1,0)[4]                    : -640.5805
 ARIMA(0,0,2)(2,1,0)[4] with drift         : -664.183
 ARIMA(0,0,2)(2,1,1)[4]                    : Inf
 ARIMA(0,0,2)(2,1,1)[4] with drift         : -661.8972
 ARIMA(0,0,3)(0,1,0)[4]                    : -642.0466
 ARIMA(0,0,3)(0,1,0)[4] with drift         : -661.4649
 ARIMA(0,0,3)(0,1,1)[4]                    : -639.9556
 ARIMA(0,0,3)(0,1,1)[4] with drift         : -664.5459
 ARIMA(0,0,3)(0,1,2)[4]                    : -638.8165
 ARIMA(0,0,3)(0,1,2)[4] with drift         : -662.5133
 ARIMA(0,0,3)(1,1,0)[4]                    : -639.977
 ARIMA(0,0,3)(1,1,0)[4] with drift         : -664.8335
 ARIMA(0,0,3)(1,1,1)[4]                    : -637.8841
 ARIMA(0,0,3)(1,1,1)[4] with drift         : -662.5329
 ARIMA(0,0,3)(2,1,0)[4]                    : -639.1253
 ARIMA(0,0,3)(2,1,0)[4] with drift         : -662.5329
 ARIMA(1,0,0)(0,1,0)[4]                    : -636.2003
 ARIMA(1,0,0)(0,1,0)[4] with drift         : -665.6442
 ARIMA(1,0,0)(0,1,1)[4]                    : -634.1766
 ARIMA(1,0,0)(0,1,1)[4] with drift         : -667.9443
 ARIMA(1,0,0)(0,1,2)[4]                    : -637.0258
 ARIMA(1,0,0)(0,1,2)[4] with drift         : -666.4038
 ARIMA(1,0,0)(1,1,0)[4]                    : -634.2753
 ARIMA(1,0,0)(1,1,0)[4] with drift         : -668.4951
 ARIMA(1,0,0)(1,1,1)[4]                    : -632.88
 ARIMA(1,0,0)(1,1,1)[4] with drift         : -666.3547
 ARIMA(1,0,0)(1,1,2)[4]                    : -650.5238
 ARIMA(1,0,0)(1,1,2)[4] with drift         : -664.1574
 ARIMA(1,0,0)(2,1,0)[4]                    : -640.0935
 ARIMA(1,0,0)(2,1,0)[4] with drift         : -666.3727
 ARIMA(1,0,0)(2,1,1)[4]                    : Inf
 ARIMA(1,0,0)(2,1,1)[4] with drift         : -664.141
 ARIMA(1,0,0)(2,1,2)[4]                    : Inf
 ARIMA(1,0,0)(2,1,2)[4] with drift         : -662.6872
 ARIMA(1,0,1)(0,1,0)[4]                    : Inf
 ARIMA(1,0,1)(0,1,0)[4] with drift         : -663.4946
 ARIMA(1,0,1)(0,1,1)[4]                    : Inf
 ARIMA(1,0,1)(0,1,1)[4] with drift         : -665.737
 ARIMA(1,0,1)(0,1,2)[4]                    : Inf
 ARIMA(1,0,1)(0,1,2)[4] with drift         : -664.1545
 ARIMA(1,0,1)(1,1,0)[4]                    : Inf
 ARIMA(1,0,1)(1,1,0)[4] with drift         : -666.2868
 ARIMA(1,0,1)(1,1,1)[4]                    : Inf
 ARIMA(1,0,1)(1,1,1)[4] with drift         : -664.1015
 ARIMA(1,0,1)(1,1,2)[4]                    : Inf
 ARIMA(1,0,1)(1,1,2)[4] with drift         : -661.8597
 ARIMA(1,0,1)(2,1,0)[4]                    : -638.1237
 ARIMA(1,0,1)(2,1,0)[4] with drift         : -664.1202
 ARIMA(1,0,1)(2,1,1)[4]                    : Inf
 ARIMA(1,0,1)(2,1,1)[4] with drift         : -661.8426
 ARIMA(1,0,2)(0,1,0)[4]                    : Inf
 ARIMA(1,0,2)(0,1,0)[4] with drift         : -662.4704
 ARIMA(1,0,2)(0,1,1)[4]                    : Inf
 ARIMA(1,0,2)(0,1,1)[4] with drift         : -669.0158
 ARIMA(1,0,2)(0,1,2)[4]                    : Inf
 ARIMA(1,0,2)(0,1,2)[4] with drift         : -666.7169
 ARIMA(1,0,2)(1,1,0)[4]                    : Inf
 ARIMA(1,0,2)(1,1,0)[4] with drift         : -668.6219
 ARIMA(1,0,2)(1,1,1)[4]                    : Inf
 ARIMA(1,0,2)(1,1,1)[4] with drift         : -666.7189
 ARIMA(1,0,2)(2,1,0)[4]                    : Inf
 ARIMA(1,0,2)(2,1,0)[4] with drift         : -666.5477
 ARIMA(1,0,3)(0,1,0)[4]                    : -647.0667
 ARIMA(1,0,3)(0,1,0)[4] with drift         : -667.3772
 ARIMA(1,0,3)(0,1,1)[4]                    : -645.7353
 ARIMA(1,0,3)(0,1,1)[4] with drift         : -669.6911
 ARIMA(1,0,3)(1,1,0)[4]                    : -646.1064
 ARIMA(1,0,3)(1,1,0)[4] with drift         : -668.2139
 ARIMA(2,0,0)(0,1,0)[4]                    : -640.7587
 ARIMA(2,0,0)(0,1,0)[4] with drift         : -663.5107
 ARIMA(2,0,0)(0,1,1)[4]                    : -640.3374
 ARIMA(2,0,0)(0,1,1)[4] with drift         : -665.7382
 ARIMA(2,0,0)(0,1,2)[4]                    : -639.9964
 ARIMA(2,0,0)(0,1,2)[4] with drift         : -664.1586
 ARIMA(2,0,0)(1,1,0)[4]                    : -640.9128
 ARIMA(2,0,0)(1,1,0)[4] with drift         : -666.287
 ARIMA(2,0,0)(1,1,1)[4]                    : -639.0462
 ARIMA(2,0,0)(1,1,1)[4] with drift         : -664.1022
 ARIMA(2,0,0)(1,1,2)[4]                    : Inf
 ARIMA(2,0,0)(1,1,2)[4] with drift         : -661.863
 ARIMA(2,0,0)(2,1,0)[4]                    : -639.8233
 ARIMA(2,0,0)(2,1,0)[4] with drift         : -664.1214
 ARIMA(2,0,0)(2,1,1)[4]                    : Inf
 ARIMA(2,0,0)(2,1,1)[4] with drift         : -661.8453
 ARIMA(2,0,1)(0,1,0)[4]                    : Inf
 ARIMA(2,0,1)(0,1,0)[4] with drift         : -662.0475
 ARIMA(2,0,1)(0,1,1)[4]                    : Inf
 ARIMA(2,0,1)(0,1,1)[4] with drift         : -667.3595
 ARIMA(2,0,1)(0,1,2)[4]                    : Inf
 ARIMA(2,0,1)(0,1,2)[4] with drift         : -665.0689
 ARIMA(2,0,1)(1,1,0)[4]                    : Inf
 ARIMA(2,0,1)(1,1,0)[4] with drift         : -667.1479
 ARIMA(2,0,1)(1,1,1)[4]                    : Inf
 ARIMA(2,0,1)(1,1,1)[4] with drift         : -667.1735
 ARIMA(2,0,1)(2,1,0)[4]                    : Inf
 ARIMA(2,0,1)(2,1,0)[4] with drift         : -667.6303
 ARIMA(2,0,2)(0,1,0)[4]                    : -659.2616
 ARIMA(2,0,2)(0,1,0)[4] with drift         : -661.5539
 ARIMA(2,0,2)(0,1,1)[4]                    : -667.5618
 ARIMA(2,0,2)(0,1,1)[4] with drift         : Inf
 ARIMA(2,0,2)(1,1,0)[4]                    : -666.0228
 ARIMA(2,0,2)(1,1,0)[4] with drift         : -668.579
 ARIMA(2,0,3)(0,1,0)[4]                    : Inf
 ARIMA(2,0,3)(0,1,0)[4] with drift         : -665.0431
 ARIMA(3,0,0)(0,1,0)[4]                    : -649.1576
 ARIMA(3,0,0)(0,1,0)[4] with drift         : -662.9122
 ARIMA(3,0,0)(0,1,1)[4]                    : -648.5549
 ARIMA(3,0,0)(0,1,1)[4] with drift         : -664.5061
 ARIMA(3,0,0)(0,1,2)[4]                    : -646.3036
 ARIMA(3,0,0)(0,1,2)[4] with drift         : -662.4044
 ARIMA(3,0,0)(1,1,0)[4]                    : -648.5489
 ARIMA(3,0,0)(1,1,0)[4] with drift         : -664.7273
 ARIMA(3,0,0)(1,1,1)[4]                    : -646.3043
 ARIMA(3,0,0)(1,1,1)[4] with drift         : -662.4341
 ARIMA(3,0,0)(2,1,0)[4]                    : -646.3013
 ARIMA(3,0,0)(2,1,0)[4] with drift         : -662.4342
 ARIMA(3,0,1)(0,1,0)[4]                    : -652.7626
 ARIMA(3,0,1)(0,1,0)[4] with drift         : -670.0218
 ARIMA(3,0,1)(0,1,1)[4]                    : Inf
 ARIMA(3,0,1)(0,1,1)[4] with drift         : -667.1187
 ARIMA(3,0,1)(1,1,0)[4]                    : Inf
 ARIMA(3,0,1)(1,1,0)[4] with drift         : -666.4217
 ARIMA(3,0,2)(0,1,0)[4]                    : Inf
 ARIMA(3,0,2)(0,1,0)[4] with drift         : -667.1087



 Best model: ARIMA(3,0,1)(0,1,0)[4] with drift         

Series: bc_ukgas_ts 
ARIMA(3,0,1)(0,1,0)[4] with drift 

Coefficients:
          ar1      ar2     ar3     ma1   drift
      -0.9933  -0.1473  0.1618  0.8785  0.0013
s.e.   0.1281   0.1397  0.1059  0.0953  0.0002

sigma^2 = 8.615e-05:  log likelihood = 341.44
AIC=-670.89   AICc=-670.02   BIC=-655.02

Note: All models will be fit to Box-Cox transformed data with \(\lambda\) = -0.445. All coefficients and terms are expressed in terms of the transformed series. Therefore, model selection, diagnostic plots, and equation reporting will all be expressed on the transformed scale. Much Later, we will back transform the model during the forecast step.

Model 1: ARIMA(3,0,1)(0,1,0)[4]

#no.constant = false because we plan to include a drift term
ukgas_model1 <-sarima(bc_ukgas_ts,3,0,1,0,1,0, 4, no.constant = FALSE) #with drift

initial  value -4.620128 
iter   2 value -4.627296
iter   3 value -4.652077
iter   4 value -4.652266
iter   5 value -4.652405
iter   6 value -4.654876
iter   7 value -4.656638
iter   8 value -4.657831
iter   9 value -4.658959
iter  10 value -4.665969
iter  11 value -4.667491
iter  12 value -4.672409
iter  13 value -4.675827
iter  14 value -4.681157
iter  15 value -4.686161
iter  16 value -4.687814
iter  17 value -4.688371
iter  18 value -4.688390
iter  19 value -4.688430
iter  20 value -4.688439
iter  21 value -4.688439
iter  22 value -4.688439
iter  22 value -4.688440
iter  22 value -4.688440
final  value -4.688440 
converged
initial  value -4.700013 
iter   2 value -4.700020
iter   3 value -4.701327
iter   4 value -4.701804
iter   5 value -4.702037
iter   6 value -4.702052
iter   7 value -4.702053
iter   7 value -4.702053
iter   7 value -4.702053
final  value -4.702053 
converged
<><><><><><><><><><><><><><>
 
Coefficients: 
         Estimate     SE t.value p.value
ar1       -0.9933 0.1281 -7.7553  0.0000
ar2       -0.1473 0.1397 -1.0544  0.2943
ar3        0.1618 0.1059  1.5279  0.1297
ma1        0.8785 0.0953  9.2151  0.0000
constant   0.0013 0.0002  5.8031  0.0000

sigma^2 estimated as 8.185962e-05 on 99 degrees of freedom 
 
AIC = -6.450844  AICc = -6.444957  BIC = -6.298283 
 

n <- length(bc_ukgas_ts)
ukgas_model1$ICs*n

      AIC      AICc       BIC 
-696.6912 -696.0554 -680.2146 

Equation of Model 1:

\[(1 + 0.9933B + 0.1473B^2 - 0.1618B^3)(1 - B^4)X_t = (1 + 0.8785B)W_t + 0.0013\] where \[W_t \sim \mathcal{N}(0,\ 8.19 \times 10^{-5})\]

We note that AR2 and AR3 coefficients, with p values of 0.2943 and 0.1297 are not significant at the .05 level.

Comments on the diagnostic plots:

There’s quite a bit of complexity to this model, with 3 AR terms, an MA term, a box-cox transformation, and a seasonal difference.

And with AR2 and AR3 terms with such high p-values, this begs the question can we achieve similar model fit with a less complex model?

One way to limit complexity is by “fixing” terms to 0. It’s kind of a guess and check procedure. We mute a term for a similar, less complex model, and do a sarima() and see if the diagnostic plots and AICc values are acceptable.

This is good practice, as it keeps the principle of parsimony in mind.

In model 2, we will observe a model that still performs well when we fix the AR2 term to 0: with good p-values on terms, diagnostic plots, and a strong AICc value

Trying less complex models by fixing coefficients to 0

Model 2: ARIMA(3,0,1)(0,1,0)[4]

Where AR2 coefficient is fixed to 0

ukgas_model2 <- sarima(bc_ukgas_ts, 3,0,1, 0,1,0, 4,
                       fixed = c(NA, 0, NA, NA, NA), #fix AR2 to 0
                       no.constant = FALSE) # With Drift

initial  value -4.620128 
iter   2 value -4.626558
iter   3 value -4.651955
iter   4 value -4.652138
iter   5 value -4.652621
iter   6 value -4.654474
iter   7 value -4.657568
iter   8 value -4.661965
iter   9 value -4.681535
iter  10 value -4.684121
iter  11 value -4.685097
iter  12 value -4.685157
iter  13 value -4.685183
iter  14 value -4.685190
iter  15 value -4.685191
iter  15 value -4.685191
iter  15 value -4.685191
final  value -4.685191 
converged
initial  value -4.695617 
iter   2 value -4.695968
iter   3 value -4.696223
iter   4 value -4.696637
iter   5 value -4.696647
iter   5 value -4.696647
iter   5 value -4.696647
final  value -4.696647 
converged
<><><><><><><><><><><><><><>
 
Coefficients: 
         Estimate     SE t.value p.value
ar1       -0.9042 0.1001 -9.0333   0e+00
ar3        0.2441 0.0704  3.4685   8e-04
ma1        0.8617 0.1050  8.2101   0e+00
constant   0.0013 0.0003  5.0138   0e+00

sigma^2 estimated as 8.277297e-05 on 100 degrees of freedom 
 
AIC = -6.459263  AICc = -6.455378  BIC = -6.332129 
 

ukgas_model2$ICs*n

      AIC      AICc       BIC 
-697.6004 -697.1808 -683.8699 

Oh WOW!!!!AICc is even lower, when i Exclude AR2, all the p values are 0 now. plots look great. THIS is my model!!! Winner winner

Equation of Model2 (with fixed AR2)

\[(1 + 0.9042B - 0.2441B^3)(1 - B^4)X_t = (1 + 0.8617B)W_t + 0.0013\] where \[W_t \sim \mathcal{N}(0,\ 8.28 \times 10^{-5})\]

Model 3: (3,0,1)(0,1,0)[4]

Where AR1 and AR2 are fixed to 0

Now let’s try to simplify further. Perhaps an AR3 term is all we need.

ukgas_model3 <- sarima(bc_ukgas_ts, 3,0,1, 0,1,0, 4,
                       fixed = c(0, 0, NA, NA, NA),
                       no.constant = FALSE)

initial  value -4.620128 
iter   2 value -4.649588
iter   3 value -4.650255
iter   4 value -4.650265
iter   5 value -4.650266
iter   5 value -4.650266
iter   5 value -4.650266
final  value -4.650266 
converged
initial  value -4.651069 
iter   2 value -4.651822
iter   3 value -4.651886
iter   4 value -4.651886
iter   4 value -4.651886
iter   4 value -4.651886
final  value -4.651886 
converged
<><><><><><><><><><><><><><>
 
Coefficients: 
         Estimate     SE t.value p.value
ar3        0.0858 0.0986  0.8696  0.3866
ma1       -0.1909 0.0961 -1.9859  0.0498
constant   0.0013 0.0002  5.8910  0.0000

sigma^2 estimated as 9.102914e-05 on 101 degrees of freedom 
 
AIC = -6.388972  AICc = -6.386665  BIC = -6.287265 
 

ukgas_model3$ICs*n

      AIC      AICc       BIC 
-690.0090 -689.7598 -679.0246 

Equation of Model 3:

\[(1 - 0.0858B^3)(1 - B^4)X_t = (1 - 0.1909B)W_t + 0.0013\] where

\[W_t \sim \mathcal{N}(0,\ 9.10 \times 10^{-5})\]

Ljung box looks poor, and 2 at lags 4 and 5 are significant. p values of coefficients deteriorated.

Therefore, we will reject this model.

Model 4: ARIMA(2,0,1)(0,1,0)[4]

ukgas_model4 <-sarima(bc_ukgas_ts, 2,0,1,0,1,0, 4, no.constant = FALSE)

initial  value -4.623576 
iter   2 value -4.625806
iter   3 value -4.648436
iter   4 value -4.648742
iter   5 value -4.648798
iter   6 value -4.649516
iter   7 value -4.650346
iter   8 value -4.651391
iter   9 value -4.652433
iter  10 value -4.652576
iter  11 value -4.652584
iter  12 value -4.652586
iter  13 value -4.652587
iter  14 value -4.652587
iter  15 value -4.652588
iter  16 value -4.652590
iter  17 value -4.652591
iter  18 value -4.652591
iter  18 value -4.652591
iter  18 value -4.652591
final  value -4.652591 
converged
initial  value -4.648542 
iter   2 value -4.648798
iter   3 value -4.649273
iter   4 value -4.649294
iter   5 value -4.649306
iter   6 value -4.649307
iter   7 value -4.649327
iter   8 value -4.649359
iter   9 value -4.649410
iter  10 value -4.649481
iter  11 value -4.649634
iter  12 value -4.649700
iter  13 value -4.649710
iter  14 value -4.649738
iter  15 value -4.649970
iter  16 value -4.650039
iter  17 value -4.650092
iter  18 value -4.650114
iter  19 value -4.650372
iter  20 value -4.650639
iter  21 value -4.651265
iter  22 value -4.651698
iter  23 value -4.652403
iter  24 value -4.652840
iter  25 value -4.652851
iter  26 value -4.652861
iter  27 value -4.652862
iter  28 value -4.652869
iter  29 value -4.652874
iter  30 value -4.652879
iter  31 value -4.652880
iter  32 value -4.652880
iter  32 value -4.652880
iter  32 value -4.652880
final  value -4.652880 
converged
<><><><><><><><><><><><><><>
 
Coefficients: 
         Estimate     SE t.value p.value
ar1        0.6776 0.1387  4.8871  0.0000
ar2        0.2231 0.0955  2.3370  0.0214
ma1       -0.8831 0.1052 -8.3927  0.0000
constant   0.0013 0.0003  4.6307  0.0000

sigma^2 estimated as 9.082442e-05 on 100 degrees of freedom 
 
AIC = -6.371728  AICc = -6.367843  BIC = -6.244594 
 

ukgas_model4$ICs*n

      AIC      AICc       BIC 
-688.1467 -687.7271 -674.4162 

Model 4 Equation:

\[(1 - 0.6776B - 0.2231B^2)(1 - B^4)X_t = (1 - 0.8831B)W_t + 0.0013\] Where \[W_t \sim \mathcal{N}(0,\ 9.08 \times 10^{-5})\]

Summary table of values (I’ve already created this in word and can send to you)

Model	Specification	AIC	AICc	All Terms Significant (p < 0.05)?	Ljung-Box Test Result
Model 1	ARIMA(3,0,1)(0,1,0)[4]	-696.691	-696.055	❌ AR2 (p = 0.29), AR3 (p = 0.13)	✅ Most p-values > 0.05 – Pass
Model 2	ARIMA(3,0,1)(0,1,0)[4], AR2 fixed to 0	-697.600	-697.181	✅ Yes	✅ Most p-values > 0.05 – Pass
Model 3	ARIMA(3,0,1)(0,1,0)[4], AR1 & AR2 fixed to 0	-690.009	-689.760	❌ AR3 (p = 0.39)	❌ Several p-values < 0.05 – Fail
Model 4	ARIMA(2,0,1)(0,1,0)[4]	-688.147	-687.727	✅ Yes	❌ Several p-values < 0.05 – Fail

So, We are going to select Model 2!!!

Forecasting using Model 2

# have to use Arima(){forecast} not arima(){stats}. 
# arima() kept dropping my drift term! Was quite frustrating
library(forecast) 

# Fit the model
ukgas.fit.model2 <- Arima(bc_ukgas_ts, 
             order = c(3,0,1), seasonal = c(0,1,0),
             fixed = c(NA, 0, NA, NA, NA),
             include.drift = TRUE)

# Check all terms are present in model, and AR2 is 0
ukgas.fit.model2$coef

         ar1          ar2          ar3          ma1        drift 
-0.904234264  0.000000000  0.244118969  0.861675905  0.001342957 

# Forecast using ARIMA model fit on Box-Cox data
ukgas.forecast.M2.bc <- forecast(ukgas.fit.model2, h = 12)

# Plot with Transformed Scale
par(mfrow = c(1,1))
plot(ukgas.forecast.M2.bc, 
     main = "ARIMA Forecast on Transformed Scale (with Fitted + History)",
     xlab = 'Year', ylab = 'BC Transformed Scale')

# Tick marks
axis(1, at = c(time(ukgas_ts)[cycle(ukgas_ts) == 1], 1986, 1987, 1988, 1989), labels = FALSE)

The plot is difficult to interpret with the Gas consumption data in this format. The data on the y axis is the Box-Cox transformed data.

A better visualization would inverse transform the y data.

# Forecast using ARIMA model fit on Box-Cox data
ukgas.forecast.M2 <- forecast(ukgas.fit.model2, h = 12)

# Inverse transform forecast components
ukgas.forecast.M2$mean   <- InvBoxCox(ukgas.forecast.M2$mean, lambda)
ukgas.forecast.M2$lower  <- InvBoxCox(ukgas.forecast.M2$lower, lambda)
ukgas.forecast.M2$upper  <- InvBoxCox(ukgas.forecast.M2$upper, lambda)
ukgas.forecast.M2$fitted <- InvBoxCox(ukgas.forecast.M2$fitted, lambda)

# Replace x (the original series used for plotting) with the true original data
ukgas.forecast.M2$x <- ukgas_ts  # ← this is key!

# Plot everything on the correct scale
par(mfrow = c(1,1))
plot(ukgas.forecast.M2, 
     main = "ARIMA Forecast on Original Scale (with Fitted + History)",
     xlab = 'Year', ylab = 'Millions of Therms (Natural Gas)')


axis(1, at = c(time(ukgas_ts)[cycle(ukgas_ts) == 1], 1986, 1987, 1988, 1989), labels = FALSE)

plot(ukgas.forecast.M2, xlim = c(1982, 1990),
     main = "ARIMA Forecast on Original Scale (with Fitted + History)",
     xlab = 'Year', ylab = 'Millions of Therms (Natural Gas)')
     
axis(1, at = c(time(ukgas_ts)[cycle(ukgas_ts) == 1], 1986, 1987, 1988, 1989), labels = TRUE)

Conclusions:

We will end here– ignore below

HOLT WINTERS Actually, this is extra. This isn’t mentioned in the prompt at all. So, lets drop it. An evaluation of 4 ARIMA models is absolutely more than enough.

HW_bc_ukgas_ts <- HoltWinters(bc_ukgas_ts) # fitting HoltWinters model to the RAW data
plot(HW_bc_ukgas_ts)

axis(1, at = time(ukgas_ts)[cycle(ukgas_ts) == 1], labels = FALSE)

forecast_HW_bc_ukgas_ts <- forecast(HW_bc_ukgas_ts)
plot(forecast_HW_bc_ukgas_ts,
     xlim = c(1982, 1990),
     ylim = c(2.00, 2.20))

axis(1, at = time(ukgas_ts)[cycle(ukgas_ts) == 1], labels = FALSE)

checkresiduals(HW_bc_ukgas_ts)

    Ljung-Box test

data:  Residuals from HoltWinters
Q* = 11.466, df = 8, p-value = 0.1767

Model df: 0.   Total lags used: 8

Since the p-value is greater than 0.05, we fail to reject the null hypothesis. This suggests there’s no significant autocorrelation left in the residuals, meaning the model fits reasonably well.