## Monday, March 21, 2011

### ARIMA, Forecasting and Python

I ported the R code found on Rob Hyndman's blog into Python + rpy2. The latter package allows calling of R code from Python which we used here to utilize the forecast package. Test R and Python codes can also be found in the zip file below.
```import rpy2.robjects as R
import rpy2.rinterface as rinterface
from rpy2.robjects.packages import importr
import rpy2.robjects.numpy2ri
from rpy2.robjects.vectors import FloatVector, StrVector
import datetime
import numpy as np
forecast = importr("forecast")
R.r.source("arima-fn.R")

n = len(x)
m = 10
terms = 4

x = np.insert(x, 10, np.nan) # put in a NaN to see what happens
n += 1

print x

R.r.assign('x',x)
R.r.assign('n',n)
R.r.assign('m',m)
R.r.assign('terms',terms)

R.r('fit = Arima(x=x, order=c(2,0,1), xreg=fourier(1:n,terms,m))')
R.r('f = forecast(fit, h=2*m, xreg=fourier(n+1:(2*m),terms,m))')
res = R.r('forecast_results(f)')
for x in res:
print x```

All Code

arima-fourier2.R

forecast_results <- function(res)
{
return (res\$mean)
}

fourier <- function(t,terms,period)
{
print (terms)
n <- length(t)
X <- matrix(,nrow=n,ncol=2*terms)
for(i in 1:terms)
{
X[,2*i-1] <- sin(2*pi*i*t/period)
X[,2*i] <- cos(2*pi*i*t/period)
}
colnames(X) <- paste(c("S","C"),rep(1:terms,rep(2,terms)),sep="")
return(X)
}

library(forecast)

y <- yy[,'price']

n <- length(y)
print (n)
m <- 10 # how many points in the future
terms = 4

fit <- Arima(y, order=c(2,0,1), xreg=fourier(1:n,terms,m))

f = forecast(fit, h=2*m, xreg=fourier(n+1:(2*m),terms,m))

print (forecast_results(f))

plot(f)

arima-fouerier.R

n <- 2000
m <- 200
y <- ts(rnorm(n) + (1:n)%%100/30, f=m)
print (y)
quit()
fourier <- function(t,terms,period)
{
n <- length(t)
X <- matrix(,nrow=n,ncol=2*terms)
for(i in 1:terms)
{
X[,2*i-1] <- sin(2*pi*i*t/period)
X[,2*i] <- cos(2*pi*i*t/period)
}
colnames(X) <- paste(c("S","C"),rep(1:terms,rep(2,terms)),sep="")
return(X)
}

library(forecast)
fit <- Arima(y, order=c(2,0,1), xreg=fourier(1:n,4,m))
f = forecast(fit, h=2*m, xreg=fourier(n+1:(2*m),4,m))
plot(f)

print (y)
print (f)

arima.py

import rpy2.robjects as R
import rpy2.rinterface as rinterface
from rpy2.robjects.packages import importr
import rpy2.robjects.numpy2ri
from rpy2.robjects.vectors import FloatVector, StrVector
import datetime
import numpy as np
forecast = importr("forecast")
R.r.source("arima-fn.R")

n = len(x)
m = 10
terms = 4

x = np.insert(x, 10, np.nan) # put in a NaN to see what happens
n += 1

print x

R.r.assign('x',x)
R.r.assign('n',n)
R.r.assign('m',m)
R.r.assign('terms',terms)

R.r('fit = Arima(x=x, order=c(2,0,1), xreg=fourier(1:n,terms,m))')
R.r('f = forecast(fit, h=2*m, xreg=fourier(n+1:(2*m),terms,m))')
res = R.r('forecast_results(f)')
for x in res:
print x

garan.dat

price
6.40000
6.33000
6.43000
6.53000
6.48000
6.33000
6.28000
6.33000
6.33000
6.18000
6.28000
6.28000
6.43000
6.48000
6.23000
6.38000
6.33000
6.23000
6.45000
6.28000
6.38000
6.43000
6.45000
6.23000
6.04000
5.84000
5.94000
5.84000
5.95000
5.89000
5.94000
6.00000
6.35000
6.13000
6.38000
6.23000
6.04000
5.79000
5.49000
5.70000