Problem with forecast se in the forecast package

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
3 messages Options
Reply | Threaded
Open this post in threaded view
|

Problem with forecast se in the forecast package

Ajay Shah
I was experimenting with the most impressive forecast package. I tried to
rig up a simplest demo of auto.arima() in action. The forecast sd looks odd
to me. Here's a demo.

## Reproducibility
set.seed(101)

## Simulate a series from AR(1)
x <- arima.sim(list(order = c(1,0,0), ar = 0.8), n = 1000)
sd(x)
  # This is 1.566 which seems right

## Using ar()
m1 <- ar(x)

## Use forecast::auto.arima()
library(forecast)
m2 <- auto.arima(x)

## How did they fare?
m1
m2
  # Both of them got the right model.

f1 <- predict(m1, n.ahead=10)
data.frame(pointest=as.numeric(f1$pred), sd=as.numeric(f1$se))
f2 <- forecast(m2, h=10)
data.frame(pointest=cbind(as.numeric(f2$mean),
sd=(f2$upper[,2]-f2$lower[,2])/2))

Let me show you the two results.

> data.frame(pointest=as.numeric(f1$pred), sd=as.numeric(f1$se))
     pointest        sd
1  -0.3186451 0.9628466
2  -0.2898804 1.2264023
3  -0.2671877 1.3649788
4  -0.2492853 1.4445294
5  -0.2351619 1.4918999
6  -0.2240198 1.5206374
7  -0.2152297 1.5382519
8  -0.2082952 1.5491136
9  -0.2028244 1.5558355
10 -0.1985085 1.5600043

This seems right. The forecast sd is small at first, as we're exploiting
the time series structure. But as we go deep into the future, the AR model
is useless and we're down to the unconditional sd which is ~ 1.566.

> f2 <- forecast(m2, h=10)
> data.frame(pointest=cbind(as.numeric(f2$mean),
sd=(f2$upper[,2]-f2$lower[,2])/2))
   pointest.V1 pointest.sd
1  -0.28090691    1.887602
2  -0.22221139    2.406792
3  -0.17578030    2.681016
4  -0.13905100    2.839174
5  -0.10999628    2.933809
6  -0.08701255    2.991505
7  -0.06883127    3.027050
8  -0.05444897    3.049081
9  -0.04307186    3.062787
10 -0.03407199    3.071333

This seems out of line. Right at n.ahead=1 the forecast sd is 1.8876 which
is > 1.566. And going beyond, the forecast sd goes up to 3.07. What am I
missing?

--
Ajay Shah
[hidden email]
http://www.mayin.org/ajayshah
http://ajayshahblog.blogspot.com

        [[alternative HTML version deleted]]

_______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-finance
-- Subscriber-posting only. If you want to post, subscribe first.
-- Also note that this is not the r-help list where general R questions should go.
Reply | Threaded
Open this post in threaded view
|

Re: Problem with forecast se in the forecast package

Adam Ginensky
Hi,
I don't know what's going on, but notice the 'bad' sd estimates are
essentially double the 'correct' sd estimates.  There is likely a
connection.

Best,

Adam Ginensky

On Mon, Dec 26, 2016 at 7:19 AM, Ajay Shah <[hidden email]> wrote:

> I was experimenting with the most impressive forecast package. I tried to
> rig up a simplest demo of auto.arima() in action. The forecast sd looks odd
> to me. Here's a demo.
>
> ## Reproducibility
> set.seed(101)
>
> ## Simulate a series from AR(1)
> x <- arima.sim(list(order = c(1,0,0), ar = 0.8), n = 1000)
> sd(x)
>   # This is 1.566 which seems right
>
> ## Using ar()
> m1 <- ar(x)
>
> ## Use forecast::auto.arima()
> library(forecast)
> m2 <- auto.arima(x)
>
> ## How did they fare?
> m1
> m2
>   # Both of them got the right model.
>
> f1 <- predict(m1, n.ahead=10)
> data.frame(pointest=as.numeric(f1$pred), sd=as.numeric(f1$se))
> f2 <- forecast(m2, h=10)
> data.frame(pointest=cbind(as.numeric(f2$mean),
> sd=(f2$upper[,2]-f2$lower[,2])/2))
>
> Let me show you the two results.
>
>> data.frame(pointest=as.numeric(f1$pred), sd=as.numeric(f1$se))
>      pointest        sd
> 1  -0.3186451 0.9628466
> 2  -0.2898804 1.2264023
> 3  -0.2671877 1.3649788
> 4  -0.2492853 1.4445294
> 5  -0.2351619 1.4918999
> 6  -0.2240198 1.5206374
> 7  -0.2152297 1.5382519
> 8  -0.2082952 1.5491136
> 9  -0.2028244 1.5558355
> 10 -0.1985085 1.5600043
>
> This seems right. The forecast sd is small at first, as we're exploiting
> the time series structure. But as we go deep into the future, the AR model
> is useless and we're down to the unconditional sd which is ~ 1.566.
>
>> f2 <- forecast(m2, h=10)
>> data.frame(pointest=cbind(as.numeric(f2$mean),
> sd=(f2$upper[,2]-f2$lower[,2])/2))
>    pointest.V1 pointest.sd
> 1  -0.28090691    1.887602
> 2  -0.22221139    2.406792
> 3  -0.17578030    2.681016
> 4  -0.13905100    2.839174
> 5  -0.10999628    2.933809
> 6  -0.08701255    2.991505
> 7  -0.06883127    3.027050
> 8  -0.05444897    3.049081
> 9  -0.04307186    3.062787
> 10 -0.03407199    3.071333
>
> This seems out of line. Right at n.ahead=1 the forecast sd is 1.8876 which
> is > 1.566. And going beyond, the forecast sd goes up to 3.07. What am I
> missing?
>
> --
> Ajay Shah
> [hidden email]
> http://www.mayin.org/ajayshah
> http://ajayshahblog.blogspot.com
>
>         [[alternative HTML version deleted]]
>
> _______________________________________________
> [hidden email] mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-finance
> -- Subscriber-posting only. If you want to post, subscribe first.
> -- Also note that this is not the r-help list where general R questions should go.

_______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-finance
-- Subscriber-posting only. If you want to post, subscribe first.
-- Also note that this is not the r-help list where general R questions should go.
Reply | Threaded
Open this post in threaded view
|

Re: Problem with forecast se in the forecast package

Ajay Shah
On 28 December 2016 at 03:20, Adam Ginensky <[hidden email]> wrote:

> Hi,
> I don't know what's going on, but notice the 'bad' sd estimates are
> essentially double the 'correct' sd estimates.  There is likely a
> connection.
>

Ouch! Got it. Thanks,

--
Ajay Shah
[hidden email]
http://www.mayin.org/ajayshah
http://ajayshahblog.blogspot.com

        [[alternative HTML version deleted]]

_______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-finance
-- Subscriber-posting only. If you want to post, subscribe first.
-- Also note that this is not the r-help list where general R questions should go.