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.

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)

yy <- read.table ("garan.dat", header=T)

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)

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") x = np.loadtxt("garan.dat", skiprows=1) 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)

yy <- read.table ("garan.dat", header=T)

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")

x = np.loadtxt("garan.dat", skiprows=1)

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