1
$\begingroup$

I have been given data to forecast however it has a negative figure within the data which then, when doing a log transformation to make the series stationary, the ARIMA script i have written won't work.

datan<-c(144627.7451,575166.2487,854245.7137,1230639.153,1160052.421,479928.7072,-261427.4238,1181746.229,168251.621,556741.5149,1840484.518,1704679.404,1878380.278,1865288.502,1849340.253,1965974.112,2093192.242,1912399.391,2633179.421,2134618.008,2070856.492,1238565.331)

freqdata<-4
startdata<-c(9,2)
horiz<-4
datats<-ts(datan,frequency=freqdata,start=startdata)
force.log<-"log"
datadates<-as.character(c("9q2","9q3","9q4","10q1","10q2","10q3","10q4","11q1","11q2","11q3","11q4","12q1","12q2","12q3","12q4","13q1","13q2","13q3","13q4","14q1","14q2","14q3"))
dataMAT<-matrix(0,ncol=freqdata,nrow=(length(datats)+freqdata),byrow=TRUE)
for (i in 1:freqdata)
  {dataMAT[,i]<-c(rep(0,length=i-1),lag(datats,k=-i+1),rep(0,length=freqdata-i+1))}
dataind<-dataMAT[c(-1:(-freqdata+1),-(length(dataMAT[,1])-freqdata+1):-(length(dataMAT[,1]))),]
dataind2<-data.frame(dataind)
lm1<-lm(X1~.,data=dataind2)
lm2<-lm(X1~X2+dataind2[,length(dataind2[1,])],data=dataind2)
library(lmtest)
library(car)
bptest1<-bptest(lm1)
bptest2<-bptest(lm2)
gqtest1<-gqtest(lm1)
ncvtest1<-ncvTest(lm1)
ncvtest2<-ncvTest(lm2)
if(force.log=="level") 
  {aslog<-"n"}else
    {{if(force.log=="log")
       {aslog<-"y"}else
         {if(bptest1$p.value<0.1|bptest2$p.value<0.1|gqtest1$p.value<0.1|ncvtest1$p<0.1|ncvtest2$p<0.1)
           {aslog<-"y"}else
              {aslog<-"n"}}}}
if(aslog=="y")
  {dataa<-log(datats)}else
    {dataa<-datats}
startLa<-startdata[1]+trunc((1/freqdata)*(length(dataa)-horiz))
startLb<-1+((1/freqdata)*(length(dataa)-horiz)-trunc((1/freqdata)*(length(dataa)-horiz)))*freqdata
startL<-c(startLa,startLb)
K<-ts(rep(dataa,length=length(dataa)-horiz),frequency=freqdata,start=startdata)
L<-ts(dataa[-1:-(length(dataa)-horiz)],frequency=freqdata,start=startL)
library(strucchange)
efp1rc<-efp(lm1,data=dataind2,type="Rec-CUSUM")
efp2rc<-efp(lm2,data=dataind2,type="Rec-CUSUM")
efp1rm<-efp(lm1,data=dataind2,type="Rec-MOSUM")
efp2rm<-efp(lm2,data=dataind2,type="Rec-MOSUM")
plot(efp2rc)
lines(efp1rc$process,col ="darkblue")
plot(efp2rm)
lines(efp1rm$process,col="darkblue")
gefp2<-gefp(lm2,data=dataind2)
plot(gefp2)
plot(dataa)
pacf(dataa)
sctest(efp2rc)
cat("log series,y/n?:",aslog)

then i want to run arima to get the forecasts

library(tseries)
library(forecast)
max.sdiff<-3
arima.force.seasonality<-"n"
kpssW<-kpss.test(dataa,null="Level")
ppW<-tryCatch({ppW<-pp.test(dataa,alternative="stationary")},error=function(ppW){ppW<-list(error="TRUE",p.value=0.99)})
adfW<-adf.test(dataa,alternative="stationary",k=trunc((length(dataa)-1)^(1/3)))
if(kpssW$p.value<0.05|ppW$p.value>0.05|adfW$p.value>0.05)
  {ndiffsW=1}else
    {ndiffsW=0}
aaW<-auto.arima(dataa,max.D=max.sdiff,d=ndiffsW,seasonal=TRUE,allowdrift=FALSE,stepwise=FALSE,trace=TRUE,seasonal.test="ch")
orderWA<-c(aaW$arma[1],aaW$arma[6],aaW$arma[2])
orderWS<-c(aaW$arma[3],aaW$arma[7],aaW$arma[4])
if(sum(aaW$arma[1:2])==0)
  {orderWA[1]<-1}else
    {NULL}
if(arima.force.seasonality=="y")
  {if(sum(aaW$arma[3:4])==0)
    {orderWS[1]<-1}else
      {NULL}}else
        {NULL}
Arimab<-Arima(dataa,order=orderWA,seasonal=list(order=orderWS),method="ML")
fArimab<-forecast(Arimab,h=8,simulate=TRUE,fan=TRUE)
if(aslog=="y")
  {fArimabF<-exp(fArimab$mean[1:horiz])}else
    {fArimabF<-fArimab$mean[1:horiz]}
plot(fArimab,main="ARIMA Forecast",sub="blue=fitted,red=actual") 
lines(dataa,col="red",lwd=2) #changes colour and size of dataa
lines(ts(append(fitted(Arimab),fArimab$mean[1]),frequency=freqdata,start=startdata),col="blue",lwd=2)
if(aslog=="y")
  {Arimab2f<-exp(fArimab$mean[1:horiz])}else
    {Arimab2f<-fArimab$mean[1:horiz]} 
start(fArimab$mean)->startARIMA
ArimaALTf<-ts(prettyNum(Arimab2f,big.interval=3L,big.mark=","),frequency=freqdata,start=startARIMA)
View(ArimaALTf,title="ARIMA2 final forecast") #brings up table of the forecasts
summary(Arimab)

If anyone can help me figure out how to forecast this data with the negative i will be really grateful!!

$\endgroup$

migrated from stackoverflow.com Sep 12 '14 at 14:53

This question came from our site for professional and enthusiast programmers.

  • $\begingroup$ Consider making a minimal reproducible example. Is all of the code you posted relevant to your question? $\endgroup$ – Will Beason Sep 12 '14 at 14:12
  • $\begingroup$ @WillBeason i put in my whole code so that people could run it and see exactly what could be changed and where as before i have been told i didnt put enough information in $\endgroup$ – Summer-Jade Gleek'away Sep 12 '14 at 14:15
1
$\begingroup$

You could shift the data by adding a constat, e.g. datats <- datats + 500000, so that all the values are positive and logs can be taken. Remember to undo this shift and recover the original level when obtaining forecasts (as you already did undoing the logarithmic transformation by taking the exponential).

Why do you take logarithms? The data do not seem to show an increasing variance. I would rather say that there is a level shift around observation $11$.

$\endgroup$
  • $\begingroup$ Thankyou for your answer! I had tried this but I wasnt sure if it would cause the forecasts to be inaccurate. I was told to force the logarithm as it created more accurate results...is this not true then? $\endgroup$ – Summer-Jade Gleek'away Sep 19 '14 at 8:33
  • $\begingroup$ I would use the original data and search for an ARIMA model including as regressor a level shift at observation 11 (a dummy with zeros at observations 1 to 10 and ones at observations 11 to the end of the sample). After the series is adjusted for the shift at time point 11 the variance is relatively homogenous, so there is no need to rescale the data. $\endgroup$ – javlacalle Sep 19 '14 at 11:22
  • $\begingroup$ How exactly would i do that? I haven't come across that before $\endgroup$ – Summer-Jade Gleek'away Sep 19 '14 at 11:25
  • 1
    $\begingroup$ Once you create the dummy, for example: d <- as.numeric(rep(0, length(datats)) + seq_along(datats)>10), you can include it through argument xreg in arima or auto.arima, e.g.: forecast::auto.arima(datats, xreg=d). You can also use the package tsoutliers, e.g. tso(datats, remove.method="bottom-up"), which detects the level shift. $\endgroup$ – javlacalle Sep 19 '14 at 12:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.