Performing a time series ARIMA model on natural gas power demand using the forecast package from R

Question

I've been attempting to forecast natural gas power demand and how it is affected by temperature and price. I'm not sure if I have done everything correctly (relatively new to R), but I do seem to get relevant data other than I can't seem to change my forecast period, nor am I sure this is an appropriate model for this data. Hopefully someone can provide me with some guidance.

Data: demand.csv

library(forecast)
data = read.csv("demand.csv")

# Create matrix of numeric predictors
xreg <- cbind(weather=data$Weather,price=data$Price,m1=data$M1,
m2=data$M2,m3=data$M3,m4=data$M4,m5=data$M5,m6=data$M6,
m7=data$M7,m8=data$M8,m9=data$M9,m10=data$M10,m11=data$M11)

# Rename columns
colnames(xreg) <- c("Weather","Price","Jan","Feb","Mar","Apr",
"May","Jun","Jul","Aug","Sep","Oct","Nov")

# Variable to be modelled
demandTS <- ts(data$Demand, frequency=12)

# Find ARIMAX model
demandArima <- auto.arima(demandTS, xreg=xreg)
demand.fcast <- forecast(demandArima, xreg=xreg)
plot(demand.fcast)

Thank you for any help.

References:

How to setup xreg argument in auto ARIMA in R From auto ARIMA to forecast in R

Answer 1

Unfortunately you have few technical errors here.

You cannot make ARIMAX-model with library(forecast) function auto.arima. Xreg argument makes it regression model with ARMA errors. That is something which I had to learn hard way by wondering the results... :)

And you have to supply FUTURE values for the xreg argument in the forecast-function. Split your data into two parts: 1) one to fit model 2) future values for the exogenous variables. Auto.arima does not forecast future values of xreg variables by ARIMA-models.

If you want to try ARIMAX-models try library(TSA) with arimax-function which of course has different syntax than auto.arima - function... :)

EDIT:

Here is example for using auto.arima with xreg argument with data set having first data for model parameters estimation and then for forecasting.

library(forecast);  
apu=read.table("demo_1.csv",sep=";",dec=",",header=TRUE);  
apux=read.table("demo_2.csv",sep=";",dec=",",header=TRUE);  
apuxx=read.table("demo_3_xreg.csv",sep=";",dec=",",header=TRUE);  
apu2=ts(data=apu[2],start=c(2011,1),deltat=1/365);  
apu3=ts(data=apu[3],start=c(2011,1),deltat=1/365);  
apu4=ts(data=apux[1],start=c(2011,1),deltat=1/365);  
Acf(apu2);  
Pacf(apu2);  
apu5=ts.intersect(apu3,apu4);  
apu6=ts(data=apuxx[3],start=c(2013,263),deltat=1/365);  
apu7=ts(data=apuxx[2],start=c(2013,263),deltat=1/365);  
apu8=ts.intersect(apu6,apu7);   


sarimax=auto.arima(apu2, d=NA, D=NA, max.p=5, max.q=5,  
     max.P=365, max.Q=365, max.order=5, start.p=2, start.q=2,  
     start.P=1, start.Q=1, stationary=FALSE, seasonal=TRUE,  
     ic=c("aicc","aic", "bic"), stepwise=TRUE, trace=FALSE,  
     approximation=(length(apu2)>100 | frequency(apu2)>12), xreg=apu5,  
     test=c("kpss","adf","pp"), seasonal.test=c("ocsb","ch"),  
     allowdrift=TRUE, lambda=0, parallel=FALSE, num.cores=NULL);  
print(sarimax$arma);  
    print(accuracy(sarimax));  
    print(sarimax$coef);  
plot(sarimax$residuals);    
    print(Box.test(sarimax$residuals,lag=30,type=c("Ljung-Box")));    

sarimaxpredicts=forecast(sarimax, h=7,level=c(75,80,90,95), fan=FALSE, xreg=apu8,  lambda=sarimax$lambda,bootstrap=FALSE, npaths=5000);  
plot(sarimaxpredicts);  

Answer 2

Unobserved Components Model (UCM) is a state space modeling approach to time series forecasting and regression analysis. It is a very flexible modeling approach, easy to interpret method. I'm not trained in statistics. UCM is a very intuitive method that a non statistician like me could easily adopt and focus on solving the problem as opposed to spending time on model building. UCM is also a very general model and every sensible ARIMA has an UCM equivalent, however not all the UCM have an equivalent ARIMA model.

I used SAS to model the UCM. Hovwever, there are other tools such as Stata, Oxmetrics that also have UCM. I don't think there is an easy to use UCM package in R. There is an excellent text An Introduction to State Space Time Series Analysis by Jacques J.F. Commandeur (Author), Siem Jan Koopman (Author) that you could refer to.

Below is the code, that I used to model, I also regressed lag value of price, becuase it is reasonable to assume that current months price can have an effect on future months demand.

data demandm;
            set demand;

            weather_1 = lag(weather); *lag 1;
            price_1 = lag(price); *lag 1;

            weather_2 = lag2(weather); *lag 2;
            price_2 = lag2(price); *lag 2;

            if weather_2 = . then delete;
            if price_2 = . then delete;
    run;

ods graphics on ; 
   proc ucm data = demandm; 
      id date interval = month; 
      model demand =weather price price_1; 
      irregular p = 1; 
      season length = 12 variance = 0 noest plot = smooth;
      estimate back=5 plot=(normal acf);
      forecast lead=5 back=5 plot=decomp; 
   run ; 
  ods graphics off ;

About the code:

• I created a dataset with lag1 and lag2 variable for weather and price and deleted the first 2 observations because it does not include 2 lags.
• I used a basic structural model in Proc UCM with only irregular component and seasonal component. In addition, I regressed demand on weather, current months price and past months price. All other variable were not significant. I also held out last 5 months to see how well this model predicts the unknown future data.
• You could forecast the data by simply changing the forecast statement within proc ucm, however you need to supply the independent variables weather and price for future values.

Below is the output:

    Irregular           Error Variance    2.416887E12     3.6646E11       6.60      <.0001
    Irregular           AR_1                  0.35428       0.10069       3.52      0.0004
    Weather             Coefficient             97713       11412.6       8.56      <.0001
    Price               Coefficient            496679      156407.4       3.18      0.0015
    price_1             Coefficient           -770025      164737.7      -4.67      <.0001

Since the analysis of residuals showed some pattern left, I modeled irregular component as an autoregressive order = 1. See the code in the irregular component where I have put p = 1. Price at lag2 and weather at both lags are not significant. With regards to significant individual variables,

• Each unit change in weather has an increase of ~98,000 units in
  demand.
• Current months price has an increase of ~496,000 units of demand.
• Current months increase in price has a ~770,000 unit decrease in next months demand.

The model also automatically identified 2 outliers in July 2012 and August 2007 which can be corrected by changing the SAS code using a dummy code variables.

Below are some nice graphical outputs (read in clockwise):

1. the most important plot is the residual auto correlation which shows that there is no pattern left.
2. The season plot where you you can see there is an increase in demand during summer and winder and low demand in spring and fall. Seems intuitive, but the great strength of UCM is that it gives you an interpretable components such as this.
3. Final plot is the model with fitted data.
4. is the out of sample forecast, where you have 5 values that have been held out and the model applied as you can visually see, the current model fits the out of sample data very well, you could use this to compare other class of models.

Hope this is helpful.

Answer 3

Fortunately one doesn't have to assume the lag structure that is appropriate as this can be suggested via the Impulse Response Weights which can be(was) modified empirically via model diagnostics. Additionally one doesn't have to assume the number of seasonal indicators and their start dates as this can be easily found via Intervention Detection. Most importantly one doesn't have to assume that the error variance is constant as the Box-Cox test can be optionally automatically evaluated for different transforms. Finally one doesn't have to manually add Pulse Indicators as modern software will perform that.

Following is a model (using scaled data ) that was automatically developed by using the freely available demo from http://www.autobox.com/30day.exe .You can download this 30 day trial version yourself and duplicate these results in a matter of minutes ( without limiting guesses as to model form ). I was one of the developers of AUTOBOX and can vouch for it's thoroughness. If you have any questions/issues contact the folks at Automatic Forecasting Systems.

The Actual/Cleansed graph illustrates the points that seasonal indicators and pulses were identified.

The table of forecasts is here

with accompanying graph

In closing the Actual/Fit and Forecasts

The first period out forecast of 80 is higher than the observed value of 53 but with a standard deviation of approximately 28 , the value of 53 is about 1 standard deviation away. Forecasters forecast fo this first period out is more accurate but as somebody once said "one swallow does not make a summer" .