# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# # Introduction

# Here is a simple approach to prediction of annual water usage where I practiced time-series forecasting.
#
# **Box-Jenkins Method is used for data analysis using non-seasonal ARIMA (non-stationary) model.**
#
# The dataset provides the annual water usage in Baltimore from 1885 to 1963, or 79 years of data.
#
# The values are in the units of liters per capita per day, and there are 79 observations.
#
# The dataset is credited to Hipel and McLeod, 1994.
#
#
# The RMSE performance measure and walk-forward validation are used for model evaluation.
#
# Literature:
#
# [Time Series Analysis and Forecasting by Example](https://books.google.si/books/about/Time_Series_Analysis_and_Forecasting_by.html?id=Bqm5kJC8hgMC&printsec=frontcover&source=kp_read_button&redir_esc=y#v=onepage&q&f=false)
#
# [How to Remove Trends and Seasonality with a Difference Transform in Python](https://machinelearningmastery.com/remove-trends-seasonality-difference-transform-python/)
#
# [Autocorrelation in Time Series Data](https://www.influxdata.com/blog/autocorrelation-in-time-series-data/)

# ## Import libraries

# +
from matplotlib import pyplot
import matplotlib.cm as cm
# %matplotlib inline

import numpy as np

from pandas import read_csv, concat, Grouper, DataFrame, datetime, Series
from pandas.plotting import lag_plot, autocorrelation_plot

import warnings
from statsmodels.tsa.arima_model import ARIMA, ARIMAResults
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from statsmodels.tsa.stattools import adfuller

from sklearn.metrics import mean_squared_error
# -

# # Data preparation

# ### Import data and split to train/test and validation set

# +
series = read_csv('water.csv', header=0, index_col=0, parse_dates=True, squeeze=True)

split_point = len(series) - 10 # how may points for validation

dataset, validation = series[0:split_point], series[split_point:]
print('Dataset: %d time points \nValidation: %d time points' % (len(dataset), len(validation)))

dataset.to_csv('dataset.csv', header=False)
validation.to_csv('validation.csv', header=False)
# -

# ### Summary statistics

# +
# summary statistics of time series
series = read_csv('dataset.csv', header=None, index_col=0, parse_dates=True, squeeze=True)
print(series.describe())

# line plot
pyplot.figure(num=0, figsize=(5*2.35,5), dpi=80, facecolor='w', edgecolor='k')
series.plot()
pyplot.xlabel('time (y)')
pyplot.ylabel('water usage (lpcd)')
pyplot.title('water usage over time')
pyplot.grid(True)
pyplot.show()

# histogram plot
pyplot.figure(num=1, figsize=(5*2,5), dpi=80, facecolor='w', edgecolor='k')

#pyplot.figure(1)
pyplot.subplot(121)
series.hist()
pyplot.xlabel('water usage (lpcd)')
pyplot.ylabel('')
pyplot.title('histogram')
pyplot.grid(True)


# density plot
pyplot.subplot(122)
series.plot(kind='kde')
pyplot.xlabel('water usage (lpcd)')
pyplot.ylabel('')
pyplot.title('density plot')
pyplot.grid(True)

pyplot.tight_layout()
pyplot.show()

# +
# boxplots of time series, the last decade is omitted

pyplot.figure(num=2, figsize=(5,5), dpi=80, facecolor='w', edgecolor='k')

series = read_csv('dataset.csv', header=None, index_col=0, parse_dates=True, squeeze=True)
groups = series.groupby(Grouper(freq='10YS'))
decades = DataFrame()

for name, group in groups:
    if len(group.values) is 10:
        decades[name.year] = group.values
        
decades.boxplot()

pyplot.xlabel('time (decade)')
pyplot.ylabel('water usage (lpcd)')
pyplot.title('boxplot, groupd by decades')
pyplot.show()
# -

# ## Persistence model - baseline RMSE

# +
# evaluate a persistence model

# load data
series = read_csv('dataset.csv', header=None, index_col=0, parse_dates=True, squeeze=True)

# prepare data
X = series.values
X = X.astype('float32')
train_size = int(len(X) * 0.50)
train, test_baseline = X[0:train_size], X[train_size:]

# walk-forward / Rolling Window / Rolling Forecast validation
history = [x for x in train]
predictions = list()

for i in range(len(test_baseline)):
    # predict
    yhat = history[-1]
    predictions.append(yhat)
    # observation
    obs = test_baseline[i]
    history.append(obs)
    #print('Predicted=%.3f, Expected=%3.f' % (yhat, obs))

# report performance
rmse = np.sqrt(mean_squared_error(test_baseline, predictions))
print('Peristence RMSE: %.3f' % rmse)

# plot
pyplot.figure(num=2, figsize=(5,5), dpi=80, facecolor='w', edgecolor='k')

pyplot.plot(test_baseline)
pyplot.plot(predictions, color='red')

pyplot.xlabel('time')
pyplot.ylabel('water usage (lpcd)')
pyplot.title('persistence model')
pyplot.grid(True)

pyplot.show()


# -

# # Manually configure ARIMA

# ## Detrend data by differencing

# +
# create and summarize a stationary version of the time series

# create a differenced series
def difference(dataset):
    diff = list()
    for i in range(1, len(dataset)):
        value = dataset[i] - dataset[i - 1]
        diff.append(value)
    return Series(diff)

series = read_csv('dataset.csv', header=None, index_col=0, parse_dates=True, squeeze=True)
X = series.values
X = X.astype('float32')

# difference data for detrending
stationary = difference(X)
stationary.index = series.index[1:]

# check if stationary
result = adfuller(stationary)
print('ADF Statistic: %f' % result[0])
print('p-value: %f' % result[1])
print('Critical Values:')
for key, value in result[4].items():
    print('\t%s: %.3f' % (key, value))
    
# plot differenced data
stationary.plot()
pyplot.title('differenced data')
pyplot.xlabel('time (y)')
pyplot.ylabel('d(water usage (lpcd)) / dt')
pyplot.show()

# save
stationary.to_csv('stationary.csv', header=False)
# -

# One step differencing (d=1) appears to be enough.

# ## Autocorrelation and partial autoorrelation
# #### estimates for lag *p* and order of MA model *q*
#
# p is the order (number of time lags) of the autoregressive model
#
# d is the degree of differencing (the number of times the data have had past values subtracted)
#
# q is the order of the moving-average model

# +
# ACF and PACF plots of the time series

series = read_csv('dataset.csv', header=None, index_col=0, parse_dates=True, squeeze=True)

pyplot.figure()
pyplot.subplot(211)
plot_acf(series, lags=20, ax=pyplot.gca())
pyplot.xlabel('lag (d)')



pyplot.subplot(212)
plot_pacf(series, lags=20, ax=pyplot.gca())
pyplot.xlabel('lag (d)')

pyplot.tight_layout()
pyplot.show()
# -

# A good starting point for the p is 4 and q as 1.

# ### Evaluate a manually configured ARIMA model

# +
# load data
series = read_csv('dataset.csv', header=None, index_col=0, parse_dates=True, squeeze=True)

# prepare dataa
X = series.values
X = X.astype('float32')
train_size = int(len(X) * 0.50)
train, test = X[0:train_size], X[train_size:]

# walk-forward validation
history = [x for x in train]
predictions = list()

for i in range(len(test)):
    # predict
    model = ARIMA(history, order=(4,1,1))
    model_fit = model.fit(disp=0)
    yhat = model_fit.forecast()[0]
    predictions.append(yhat)
    # observation
    obs = test[i]
    history.append(obs)
    #print('>Predicted=%.3f, Expected=%3.f' % (yhat, obs))
    
# report performance
rmse = np.sqrt(mean_squared_error(test, predictions))
print('RMSE: %.3f' % rmse)


# -

# Worse performance than baseline (persistence) model!

# # Grid search ARIMA parameters

# ### Define functions

# +
# evaluate an ARIMA model for a given order (p,d,q) and return RMSE
def evaluate_arima_model(X, arima_order):
    
    # prepare training dataset
    X = X.astype('float32')
    train_size = int(len(X) * 0.50)
    train, test = X[0:train_size], X[train_size:]
    history = [x for x in train]
    
    # make predictions
    predictions = list()
    for t in range(len(test)):
        model = ARIMA(history, order=arima_order)
        
        # model_fit = model.fit(disp=0)
        model_fit = model.fit(trend='nc', disp=0) # disable the automatic addition of a trend constant
        yhat = model_fit.forecast()[0]
        predictions.append(yhat)
        history.append(test[t])
        
    # calculate out of sample error
    rmse = np.sqrt(mean_squared_error(test, predictions))
    return rmse

# evaluate combinations of p, d and q values for an ARIMA model
def evaluate_models(dataset, p_values, d_values, q_values):
    
    dataset = dataset.astype('float32')
    best_score, best_order = float("inf"), None
    
    for p in p_values:
        for d in d_values:
            for q in q_values:
                
                order = (p,d,q)
                
                try:
                    rmse = evaluate_arima_model(dataset, order)
                    
                    if rmse < best_score:
                        best_score, best_order = rmse, order
                        
                    print('ARIMA%s - RMSE = %.3f' % (order, rmse))
                    
                except:
                    continue
                    
    print('\nBest: ARIMA%s - RMSE = %.3f' % (best_order, best_score))
    
    return best_order


# -

# ### Run on dataset

def grid_search(series):

    # evaluate parameters
    p_values = range(0, 3)
    d_values = range(0, 2)
    q_values = range(0, 3)

    warnings.filterwarnings("ignore")

    best_order = evaluate_models(series.values, p_values, d_values, q_values)
    
    return best_order


# +
# load dataset
series = read_csv('dataset.csv', header=None, index_col=0, parse_dates=True, squeeze=True)

# search
best_order = grid_search(series)


# -

# ### Summarize residual errors - *walk-forward validation*

def residual_stats(series, best_order, bias=0):
    
    print('-----------------------------')
    
    # prepare data
    X = series.values
    X = X.astype('float32')
    train_size = int(len(X) * 0.50)
    train, test = X[0:train_size], X[train_size:]

    # walk-forward validation
    history = [x for x in train]
    predictions = list()

    for i in range(len(test)):

        # predict
        model = ARIMA(history, order=best_order)
        model_fit = model.fit(trend='nc', disp=0)
        yhat = bias + float(model_fit.forecast()[0])
        predictions.append(yhat)

        # observation
        obs = test[i]
        history.append(obs)

    #report performance
    rmse = np.sqrt(mean_squared_error(test, predictions))
    print('RMSE: %.3f\n' % rmse)

    # errors
    residuals = [test[i]-predictions[i] for i in range(len(test))]
    residuals = DataFrame(residuals)
    residuals_mean = residuals.mean()

    print('RESIDUAL STATISTICS \n')
    print(residuals.describe())

    # plot
    pyplot.figure(num=None, figsize=(5,5), dpi=80, facecolor='w', edgecolor='k')

    pyplot.subplot(211)
    residuals.hist(ax=pyplot.gca())
    pyplot.xlabel('residual')
    pyplot.ylabel('')
    pyplot.title('histogram')
    pyplot.grid(True)

    pyplot.subplot(212)
    residuals.plot(kind='kde', ax=pyplot.gca())
    pyplot.xlabel('residual')
    pyplot.ylabel('')
    pyplot.title('density plot')
    pyplot.grid(True)

    pyplot.tight_layout()
    pyplot.show()
    
    return residuals_mean


# +
series = read_csv('dataset.csv', header=None, index_col=0, parse_dates=True, squeeze=True)

residuals_mean = residual_stats(series, best_order)
# -

# We can see that the
# distribution has a right shift and that the mean is non-zero at around 1.0. This is perhaps a sign
# that the predictions are biased.
#
# The distribution of residual errors is also plotted. The graphs suggest a Gaussian-like
# distribution with a longer right tail, providing further evidence that perhaps a power transform
# might be worth exploring.
#
# We could use this information to bias-correct predictions by adding the mean residual error
# of 1.081624 to each forecast made.

# ### Make bias corrected forecasts
#

_ = residual_stats(series, best_order, bias = residuals_mean)


# Not much of an improvement after bias correction.

# # Save finalized model to file

# +
# monkey patch around bug in ARIMA class
def __getnewargs__(self):
    return ((self.endog),(self.k_lags, self.k_diff, self.k_ma))

ARIMA.__getnewargs__ = __getnewargs__

def save_model(best_order, model_name):

    # load data
    series = read_csv('dataset.csv', header=None, index_col=0, parse_dates=True, squeeze=True)

    # prepare data
    X = series.values
    X = X.astype('float32')

    # fit model
    model = ARIMA(X, order=best_order)
    model_fit = model.fit(trend='nc', disp=0)

    # bias constant
    bias = residuals_mean

    # save model
    model_fit.save(f'model_{model_name}.pkl')
    np.save(f'model_bias_{model_name}.npy', [bias])


# -

save_model(best_order_box, model_name='simple')


# # Load and evaluate the finalized model on the validation dataset

def validate_models(model_names):
    
    # load train dataset
    dataset = read_csv('dataset.csv', header=None, index_col=0, parse_dates=True, squeeze=True)
    X = dataset.values.astype('float32')
    history = [x for x in X]

    # load validation dataset
    validation = read_csv('validation.csv', header=None, index_col=0, parse_dates=True, squeeze=True)
    y = validation.values.astype('float32')
    
    # plot
    pyplot.figure(num=None, figsize=(5,5), dpi=80, facecolor='w', edgecolor='k')
    pyplot.plot(y, color=cm.Set1(0), label='actual')
    
    
    for ind, model_name in enumerate(model_names):
        
        # load model
        model_fit = ARIMAResults.load(f'model_{model_name}.pkl')
        bias = np.load(f'model_bias_{model_name}.npy')

        # make first prediction
        predictions = np.ones(len(y))
        yhat = bias + float(model_fit.forecast()[0])
        predictions[0] = yhat
        history.append(y[0])
        #print('>Predicted=%.3f, Expected=%3.f' % (yhat, y[0]))

        # rolling forecasts
        for i in range(1, len(y)):

            # predict
            model = ARIMA(history, order=(2,1,0))
            model_fit = model.fit(trend='nc', disp=0)
            yhat = bias + float(model_fit.forecast()[0])
            predictions[i] = yhat

            # observation
            obs = y[i]
            history.append(obs)
        #    print('>Predicted=%.3f, Expected=%3.f' % (yhat, obs))

        rmse = np.sqrt(mean_squared_error(y, predictions))
        print(f'RMSE {model_name}: %.3f' % rmse)

        pyplot.plot(predictions, color=cm.Set1(ind+1), label=f'{model_name} predict')

    pyplot.xlabel('time (d)')
    pyplot.ylabel('water usage (lpcd)')
    pyplot.title('Validation')

    pyplot.legend()
    pyplot.grid(True)

    pyplot.show()

validate_models(model_names=['simple'])

# # Comparison of detrend approaches

# ### Linear detrend & Box-Cox transform

# +
from statsmodels.tsa.tsatools import detrend
from scipy.stats import boxcox

figsize = (8,4)

series = read_csv('dataset.csv', header=0, index_col=0, parse_dates=True, squeeze=True)
#print(series.describe())
pyplot.figure(num=None, figsize=figsize, dpi=80, facecolor='w', edgecolor='k')

pyplot.subplot(121)
series.plot(color='k', label='data')
pyplot.subplot(122)
series.plot(kind='kde', color='k', label='data')

pyplot.xlabel('')
pyplot.ylabel('density')
pyplot.title('density plot')

pyplot.legend()
pyplot.show()


# --- linear detrend ---
series_linear = detrend(series)
#print(result.describe())
pyplot.figure(num=None, figsize=figsize, dpi=80, facecolor='w', edgecolor='k')

pyplot.subplot(121)
series_linear.plot(color='k', label='linear detrend')
pyplot.subplot(122)
series_linear.plot(kind='kde', color='k', label='linear detrend')

pyplot.xlabel('')
pyplot.ylabel('density')
pyplot.title('density plot')

pyplot.legend()
pyplot.show()

# --- Box-Cox transform ---
series_boxcox, lam = boxcox(series)
series_boxcox = Series(series_new)

# plot
pyplot.figure(num=None, figsize=figsize, dpi=80, facecolor='w', edgecolor='k')

pyplot.subplot(121)
series_boxcox.plot(color='k', label='Box-Cox')
pyplot.subplot(122)
series_boxcox.plot(kind='kde', color='k', label='Box-Cox')

pyplot.xlabel('')
pyplot.ylabel('density')
pyplot.title('density plot')

pyplot.legend()
pyplot.show()
# -

best_order_simple = grid_search(series)
best_order_linear = grid_search(series_linear)
best_order_boxcox = grid_search(series_boxcox)

_ = residual_stats(series, best_order)
_ = residual_stats(series_linear, best_order)
_ = residual_stats(series_boxcox, best_order)

# +
save_model(best_order_simple, model_name='simple')
save_model(best_order_linear, model_name='linear')
save_model(best_order_boxcox, model_name='boxcox')

validate_models(model_names = ['simple', 'linear', 'boxcox'])
# -

# Predictions with linear detrend and Box-Cox transform have lower RMSE, probably statistically insignificant.

StarcoderdataJupyterScriptsDedupFiltered