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

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)

# 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:

• you could simplify the model by excluding seasonal dummies, auto arima automatically fits seasonal models, unless you want to determine the effect of each month. Are you open to other models ? if so you could try unobserved components model. – forecaster Feb 11 '14 at 23:25
• @forecaster I am open to other models. Could you suggest some good links to information about the unobserved components model? – UnNatural Feb 12 '14 at 17:09

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);
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,
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);

• Are you sure that you cannot create an ARIMAX model using auto.arima? I am quite sure that you can. – Matteo De Felice Feb 12 '14 at 15:46
• @Analyst I'm not exactly certain how to split my data data into two parts using xreg. Could you possibly give me an example or point me to a link? – UnNatural Feb 12 '14 at 17:11
• @MatteoDeFelice here is link to the discussion by author of forecast package robjhyndman.com/hyndsight/arimax – Analyst Feb 13 '14 at 5:31
• For me ARIMAX is same as transfer function model and not regression model with ARMA structure in residuals. – Analyst Feb 13 '14 at 5:32
• I don't really completely understand the point. I assumed that an ARIMAX would be an ARIMA plus the 'xreg' part as here: stackoverflow.com/questions/15681529/… – Matteo De Felice Feb 14 '14 at 15:52

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);
run ;
ods graphics off ;


• 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.

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" .

• I appreciate the response, but I work with R due to money limitations. Thanks. – UnNatural Feb 24 '14 at 18:31
• Pehaps my results might provide a beacon for your work showing the kind of model that is possible. regards . – IrishStat Feb 25 '14 at 2:47