Python Statsmodel Arima Start [stationarity]

I just began working on time series analysis using statsmodels. I have a dataset with dates and values (for about 3 months). I am facing some issues with providing the right order

Solution 1:

To induce stationarity:

1. de-seasonalize (remove seasonality)
2. de-trend (remove trend)

There are several ways to achieve stationarity of a time series - Box-Cox family of transformations, Differencing etc., Choice of method depends on the data. Below are the commonly used tests for stationarity.

Tests for stationarity: 1. Augmented Dickey-Fuller test 2. KPSS test KPSS python code

Solution 2:

You can use R script instead statmodels. R is more powerful for statistical estimation.

If you want use python, you can run R-script from a python through os interface:

for example R script for arima estimation "arimaestimation.r":

library(rjson)

args <- commandArgs(trailingOnly=TRUE)

jsonstring = ''for(i in seq(0, length(args))) {
    if ( length(args[i]) && args[i]=='--jsondata' ) {
        jsonstring = args[i+1]
    }
}

jsonobject = fromJSON(jsonstring)
data = as.numeric(unlist(jsonobject['data']))
p = as.numeric(unlist(jsonobject['p']))
d = as.numeric(unlist(jsonobject['d']))
q = as.numeric(unlist(jsonobject['q']))

estimate = arima(data, order=c(p, d, q))

phi = c()
if (p>0) {
    for (i in seq(1, p)) {
        phi = c(phi, as.numeric(unlist(estimate$coef[i])))
    }
}
theta = c()
if (p+1 <= p+q) {
    for (i in seq(p+1, p+q)) {
        theta = c(theta, as.numeric(unlist(estimate$coef[i])))
    }
}
if (d==0) {
    intercept = as.numeric(unlist(estimate$coef[p+q+1]))
} else {
    intercept = 0.0
}

if (length(phi)) {
    if (length(phi)==1) {
        phi = list(phi)
    }
} else {
    phi = list()
}

if (length(theta)) {
    if (length(theta)==1) {
        theta = list(-1 * theta)
    } else {
        theta = -1 * theta
    }
} else {
    theta = list()
}

arimapredict = predict(estimate, n.ahead = 12)
prediction = as.numeric(unlist(arimapredict$pred))
predictionse = as.numeric(unlist(arimapredict$se))

response = list(phi=phi,
                theta=theta,
                intercept=intercept,
                sigma2=estimate$sigma2,
                aic=estimate$aic,
                prediction=prediction,
                predictionse=predictionse)

cat(toJSON(response))

And call him with python through json interface:

Rscript arima/arimaestimate.r --jsondata '{"q": 2, "p": 2, "data": [247.0, 249.0, 213.0, 154.0, 122.0, 164.0, 141.0, 174.0, 281.0, 141.0, 159.0, 168.0, 243.0, 261.0, 211.0, 303.0, 308.0, 239.0, 237.0, 185.0], "d": 1}'

and you get the answer:

{
    "phi": [],
    "theta": [
        0.407851844478153
    ],
    "intercept":0,
    "sigma2":3068.29837379914,
    "aic":210.650287294343,
    "prediction": [
        210.184175597721,
        210.184175597721,
        210.184175597721,
        210.184175597721,
        210.184175597721,
        210.184175597721,
        210.184175597721,
        210.184175597721,
        210.184175597721,
        210.184175597721,
        210.184175597721,
        210.184175597721
    ]
}