[R] predict and arima

[bryan]

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)
> > }
> 
> ...
> 
>