Indeed, there was a bug ... my current play code looks like this ...
get.best.arima <- function(x.ts, maxord=c(3,3,3,3))
{
# function based on 'Introductory Time Series with R'
# ... try and fit the best ARIMA(p,d,q,P,D,Q) model
# using all permutations from 0 to maxord for p,q,P,Q
# Assume D=1 and select d using ndiffs and the KPSS test
# for stationarity.
# ... if no model can be found - return NULL
# - let's assume D is seasonal ... D <- 1
D <- 1
# let's test for stationarity ... to determine d
pass <- FALSE
tryCatch(
{
d <- ndiffs(x.ts, alpha=0.09, test="kpss")
pass <- TRUE
}, interrupt = function(ex)
{
cat("An interrupt was detected.\n");
print(ex);
}, error = function(ex)
{
cat("An error was detected.\n");
print(ex);
}) #tryCatch
if(!pass) return(NULL)
# let's iterate over the possibilities for p, q, P, and Q
best.aic <- Inf # a big number
n <- length(x.ts)
found <- FALSE
for(p in 0:maxord[1]) for(q in 0:maxord[2])
for(P in 0:maxord[3]) for(Q in 0:maxord[4])
{
tryCatch(
{
fit <- arima(x.ts, order=c(p,q,d),
seasonal=list(order=c(P,D,Q),
period=frequency(x.ts)),
method='CSS')
fit.aic <- -2 * fit$loglik + (log(n) +
1) * length(fit$coef)
if(fit.aic < best.aic)
{
best.aic <- fit.aic
best.fit <- fit
best.model <- c(p,d,q,P,D,Q)
found <- TRUE
}
}, interrupt = function(ex)
{
cat("An interrupt was detected.\n");
print(ex);
}, error = function(ex)
{
cat("An error was detected.\n");
print(ex);
}) #tryCatch
}
if(!found) return(NULL) # ouch!
# package up and return it in a named list ...
return(list(model=best.model, fit=best.fit, aic=best.aic))
}
On Tue, 2011-02-08 at 18:40 +0100, Jose-Marcio Martins da Cruz wrote:
> I think ther's a bug here :
>
> bdp wrote:
> >
> > Some code I have been playing with to do this follows ...
> >
> > get.best.arima<- function(x.ts, minord=c(0,0,0,0,0,0),
> > maxord=c(2,1,1,2,1,1))
> > {
> > # function based on 'Introductory Time Series with R'
> > best.aic<- 1e8 # a big number
> > n<- length(x.ts)
> > for(p in minord[1]:maxord[1]) for(d in minord[2]:maxord[2]) for(q in
> > minord[3]:maxord[3])
> > {
> > for(P in minord[4]:maxord[4]) for(D in minord[5]:maxord[5])
> > for(Q in
> > minord[6]:maxord[6])
> > {
> > fit<- arima(x.ts, order=c(p,q,d),
> > seas=list(order=c(P,D,Q),
>
> maybe it should be :
>
> fit<- arima(x.ts, order=c(p,d,q),
> seas=list(order=c(P,D,Q),
>
> exchange q and d !
>
> > frequency(x.ts)), method='CSS')
> > fit.aic<- -2 * fit$loglik + (log(n) + 1) *
> > length(fit$coef)
> > if(fit.aic< best.aic) # probably should do other
> > tests here before
> > accepting
> > {
> > best.aic<- fit.aic
> > best.fit<- fit
> > best.model<- c(p,d,q,P,D,Q)
> > }
> > }
> > }
> > #print(best.aic)
> > #print(best.model)
> > return(best.fit)
> > }
>
> ...
>
>
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
