Alternative for r's Exponential smoothing state space model in python/scikit/numpy

Question

In R we have one good forecasting model like:

ets(y, model="ZZZ", damped=NULL, alpha=NULL, beta=NULL, gamma=NULL, 

phi=NULL, additive.only=FALSE, lambda=NULL, 

lower=c(rep(0.0001,3), 0.8), upper=c(rep(0.9999,3),0.98), 

opt.crit=c("lik","amse","mse","sigma","mae"), nmse=3, 

bounds=c("both","usual","admissible"), ic=c("aicc","aic","bic"),

restrict=TRUE, allow.multiplicative.trend=FALSE, use.initial.values=FALSE, ...)

In this method if we assign any variable, it automatically gets season type,trend & error type like model="ZZZ"/"AMA"/"MMZ" and some of the factors are auto adjusted to get accurate results.

• In python do we have anything similar to ets in either pandas/numpy/scipy/scikit?

  By my research:
  Ewma in pandas is similar, but we need to hardcode all the parameters to fixed ones.
  In Holtwinter we need to write detailed methods for all the trend and season types.

• So instead of that do we have any ready-made functions which takes dataframes as input and provides forecasting values, without writing any inner functions for parameters ourselves?

• Any Fine tuned regression models scikit/statsmodels?

Answer 1

After searching around a bit, I haven't found anything that seems really promising as an ets alternative for python. There are some tries though: StatsModels and pycast's Forecasting methods, which you can check if they suit your needs.

One option you can use to workaround the missing implementation is to run an R script from python using the subprocess module. There is a very good article on how to do this here.

In order to do the later:

1. You need to create an R script (ex. my_forecast.R), which will calculate (using ets) and print the forecasts on a file, or on stdout (with the cat() command), in order to use them after the script runs.

2. You can run the R script from a python script as follows:

   import subprocess
   
   # You need to define the command that will run the Rscript from the subprocess
   command = 'Rscript'
   path2script = 'path/to/my_forecast.R'
   cmd = [command, path2script]
   
   # Option 1: If your script prints to a file
   subprocess.run(cmd)
   f = open('path/to/created/file', 'r')
   (...Do stuff from here...)
   
   # Option 2: If your script prints to stdout
   forecasts = subprocess.check_output(cmd, universal_newlines=True)
   (...Do stuff from here...)
   

   You can also add arguments to your cmd, which will be utilized by your Rscript as command line arguments, as follows:

   args = [arg0, arg1, ...]
   
   cmd = [command, path2script] + args 
   Then pass cmd to the subprocess
   

---

EDIT:

I found an exemplary series of articles on the Holt-Winters Forecasting: part1, part2 and part3. Besides the easy to comprehend analysis in those articles, Gregory Trubetskoy (the author) has provided the code he developed:

Initial trend:

def initial_trend(series, slen):
    sum = 0.0
    for i in range(slen):
        sum += float(series[i+slen] - series[i]) / slen
    return sum / slen

# >>> initial_trend(series, 12)
# -0.7847222222222222

Initial Seasonal Components:

def initial_seasonal_components(series, slen):
    seasonals = {}
    season_averages = []
    n_seasons = int(len(series)/slen)
    # compute season averages
    for j in range(n_seasons):
        season_averages.append(sum(series[slen*j:slen*j+slen])/float(slen))
    # compute initial values
    for i in range(slen):
        sum_of_vals_over_avg = 0.0
        for j in range(n_seasons):
            sum_of_vals_over_avg += series[slen*j+i]-season_averages[j]
        seasonals[i] = sum_of_vals_over_avg/n_seasons
    return seasonals

# >>> initial_seasonal_components(series, 12)
# {0: -7.4305555555555545, 1: -15.097222222222221, 2: -7.263888888888888,
#  3: -5.097222222222222,  4: 3.402777777777778,   5: 8.069444444444445,  
#  6: 16.569444444444446,  7: 9.736111111111112,   8: -0.7638888888888887,
#  9: 1.902777777777778,  10: -3.263888888888889, 11: -0.7638888888888887}

Finally the algorithm:

def triple_exponential_smoothing(series, slen, alpha, beta, gamma, n_preds):
    result = []
    seasonals = initial_seasonal_components(series, slen)
    for i in range(len(series)+n_preds):
        if i == 0: # initial values
            smooth = series[0]
            trend = initial_trend(series, slen)
            result.append(series[0])
            continue
        if i >= len(series): # we are forecasting
            m = i - len(series) + 1
            result.append((smooth + m*trend) + seasonals[i%slen])
        else:
            val = series[i]
            last_smooth, smooth = smooth, alpha*(val-seasonals[i%slen]) + (1-alpha)*(smooth+trend)
            trend = beta * (smooth-last_smooth) + (1-beta)*trend
            seasonals[i%slen] = gamma*(val-smooth) + (1-gamma)*seasonals[i%slen]
            result.append(smooth+trend+seasonals[i%slen])
    return result

# # forecast 24 points (i.e. two seasons)
# >>> triple_exponential_smoothing(series, 12, 0.716, 0.029, 0.993, 24)
# [30, 20.34449316666667, 28.410051892109554, 30.438122252647577, 39.466817731253066, ...

You can place those in a file, ex: holtwinters.py inside a folder with the following structure:

forecast_folder
|
└── __init__.py
|
└── holtwinters.py

From here on, this is a python module which you can place inside every project structure you want and use it anywhere inside that project, simply by importing it.

Answer 2

ETS is now available in python as part of the Statsmodels' most recent release. Here is a simple example of its use:

from statsmodels.tsa.api import ExponentialSmoothing

# make our own periodic data
# 3 weeks of hourly data
m = 24
days = 21
x = np.array(range(days * m)) / (m / 2 / np.pi)
st = 4 * np.sin(x)
lt = np.linspace(0, 15, days * m)
bt = 100
e = np.random.normal(scale=0.5, size=x.shape)

# make the ts we wish to forecast
y = lt + st + bt + e

# our guessed parameters
alpha = 0.4
beta = 0.2
gamma = 0.01

# initialise model
ets_model = ExponentialSmoothing(y, trend='add', seasonal='add', 
seasonal_periods=24)
ets_fit = ets_model.fit(smoothing_level=alpha, smoothing_slope=beta,
smoothing_seasonal=gamma)

# forecast p hours ahead
p_ahead = 48
yh = ets_fit.forecast(p_ahead)

# plot the y, y_smoothed and y_hat ts'
plt.plot(y, label='y')
plt.plot(ets_fit.fittedvalues, label='y_smooth')
plt.plot(range(days * m, days * m + p_ahead),yh, label='y_hat')

plt.legend()
plt.show()

Here are more examples in notebook form.

Finally, here is the source code, should you wish to take a look over it.