Motivation A time-series represents the sequence of a variable recorded with time over regular intervals of time. There are various applications of time series like Stock Market Prediction (with daily and hourly data), Annual Sales/Revenue, Seasonal Temperatures changes, monthly Cloud Infrastructure cost and second level prediction of Python has 2 libraries StatsModels and Pyramid that helps to build forecasting models and predict values at a future time. As both these libraries provide forecasts for similar models namely, AR (Auto-Regressive)MA (Moving Average)ARIMA (Auto-Regressive Integrated Moving Average)SARIMA (Seasonal Auto-Regressive Integrated Moving Average) it will be interesting to see the features provided by each of the libraries along with the pros and cons of using each of them. In this blog, I try to summarise the functionalities of both of these libraries by demonstrating the Number of Active Cases for Covid-19 for any Indian state. We keep our scope limited to univariate time series analysis. DataSet The scope of the data is fro 8th March-1st May where data from 8th March-14th April has been used for training and 15th March-1st May has been used for validation and prediction. Here, we have demonstrated using state Maharashtra, however all states data and modeling can be done using the processed data from . https://github.com/sharmi1206/covid-19-analysis STATSMODEL It primarily comprises of three main modules: – Comprises of functions useful for univariate time series analysis – autoregressive models (AR), vector autoregressive models (VAR), and univariate autoregressive moving average models (ARMA), along with Non-linear models like dynamic regression and autoregression. In addition, it also provides autocorrelation, partial autocorrelation function, periodogram, and other statistical properties related to ARMA, and methods that process moving average lag-polynomials. tsa –  Equipped with functionality to model observations at any time instant with unobserved state vectors and irregular components. It has added flexibility to allow estimation with missing observations, forecasting, impulse response functions, etc. This flexibility helps it to differentiate from tsa by estimating additive and multiplicative seasonal effects on models (both on univariate and multi-variate), as well as arbitrary trend polynomials. statespace -Allows simultaneous modeling and analyzing multiple time series. vector_ar ARIMA MODELS An ARIMA model is composed of 3 constituent units which are : : . This part explores any dependent relationship between an observation and some number of lagged variables. AR Autoregression : . This part aims to make time-series stationery by subtracting or differencing an observation from observation at the previous time step of the same time series. I Integrated : This part explores the relationship between an observation and a residual error by application of moving average to lagged observations, with any given time window. MA Moving Average . In order to design a model (say ARIMA) by applying all of the above stated procedures we need to enable all parameters. We can also enable partial parameters and select to choose either AR or MA model. The parameters can be summarized as follows: : The number of lag observations included in the model, also called the lag order. p : The number of times that the raw observations are differenced also called the degree of difference. d : The size or the order of the moving average window. q To start with and in order to have different flavours of ARIMA model, the library imports are: statsmodels.tsa.ar_model AR statsmodels.tsa.arima_model ARIMA statsmodels.tsa.statespace.sarimax SARIMAX statsmodels.graphics.tsaplots plot_acf,plot_pacf from import from import from import from import MANUAL MODEL SELECTION ARIMA The following diagram illustrates the mechanism of using AR/MA/ARIMA/SARIMA by self choosing the model orders (p, q, d) and self tuning it. AR Model This Auto-Regressive model estimates a dependent variable at any timing instant based on some number of lagged observations. The variables can be summarized as follows and formulated by the below-given equation: Yt Y(t-1): Lag 1 of the time series β: Coefficient of lag1 α: Intercept term Dependence of Lag in AR Estimation, where β and α are model estimates Number of Active cases for Covid-19 can be formulated using AR model as follows. Here we put p, i.e. number of lag observations as 2, while d and q are set to 0. model_ar = ARIMA(trainActiveCases, order=( , , ))  
model_ar_fit = model_ar.fit()  

prediction_ar=model_ar_fit.forecast(len(validActiveCases))[ ]

print(model_ar_fit.summary()) residuals = pd.DataFrame(model_ar_fit.resid)
residuals.plot()
plt.show()
residuals.plot(kind= )
plt.show()
print(residuals.describe()) 2 0 0 0 # plot residual errors 'kde' The summary and the predicted output obtained from the model are illustrated below: AR Moel Summary Predicted vs Actual AR MA Model The MA model also estimates a dependent variable Yt based on the lagged forecast errors. MA Estimation with Lagged Errors where the error terms are the errors introduced from respective lags. The errors Yt and Y(t-1) are the errors from the following equations : Number of Active cases for Covid-19 can be formulated using MA model as follows. Here we put ma (moving average window), q as 2, while d and p are set to 0. model_ma = ARIMA(trainActiveCases, order=( , , ))  
model_ma_fit = model_ma.fit()  

prediction_ma=model_ma_fit.forecast(len(validActiveCases))[ ]

print(model_ma_fit.summary()) residuals = pd.DataFrame(model_ma_fit.resid)
residuals.plot()
plt.show()
residuals.plot(kind= )
plt.show()
print(residuals.describe())

model_scores.append(np.sqrt(mean_squared_error(validActiveCases,prediction_ma)))
print( ,np.sqrt(mean_squared_error(validActiveCases,prediction_ma))) 0 0 2 0 # plot residual errors 'kde' "Root Mean Square Error for MA Model: " The summary and the predicted output obtained from the model are illustrated below: MA Model Summary Predicted vs Actual MA ARIMA Model An ARIMA model comes along with both the flavours of AR and MA model where the time series is differenced at least once to make it stationary. The equation can be formulated as follows, where we see the term β coming from AR and ε–the error terms coming from MA model . ARIMA – Differencing + Combination of Parameters from AR and MA model Number of Active cases for Covid-19 can be formulated using MA model as follows. Here we put q i.e the ma (moving average window), d (order of differencing ) and p (number of lag observations) are set to 1. The code snippet for this is : model_arima = ARIMA(trainActiveCases, order=( , , ))  
model_arima_fit = model_arima.fit()  


prediction_arima=model_arima_fit.forecast(len(validActiveCases))[ ]
valid_ml[ ]=list(np.exp(prediction_arima))

print(model_arima_fit.summary()) residuals = pd.DataFrame(model_arima_fit.resid)
residuals.plot()
plt.show()
residuals.plot(kind= )
plt.show()
print(residuals.describe() 1 1 1 0 "ARIMA Model Prediction" # plot residual errors 'kde' Predicted vs Actual ARIMA We can see it’s better to use ARIMA than AR or MA model, as this is the model that closely resembles the actual outcome. SARIMAX Model SARIMAX model is built by extending the ARIMA model, discussed in our previous section. In addition to terms and terms, there are four seasonal elements that are not part of ARIMA that must be configured; they are: AR, I MA : Seasonal autoregressive order. P : Seasonal difference order. D : Seasonal moving average order. Q : The number of time steps for a single seasonal period. m Number of Active cases for Covid-19 can be formulated using the model as follows. Here we put ma (moving average window), d (order of differencing ) and p (number of lag observations) are set to 1. We also enforce seasonal order as (0,0,0,12) to study its behavior, even though Covid-19 does not show any seasonality, evident from research and daily reports. SARIMAX The code snippet for SARIMAX (without exogenous variable) can be given us: model_sarima = sm.tsa.statespace.SARIMAX(trainActiveCases, order=( , , ), seasonal_order=( , , , ),
                                 enforce_stationarity= , enforce_invertibility= )
model_sarima_fit = model_sarima.fit()  


prediction_sarima=model_sarima_fit.forecast(np.shape(validActiveCases)[ ])
valid_ml[ ]=list(np.exp(prediction_sarima))

print(model_sarima_fit.summary()) residuals = pd.DataFrame(model_sarima_fit.resid)
residuals.plot()
plt.show()
residuals.plot(kind= )
plt.show()
print(residuals.describe()) 1 1 1 0 0 0 0 False False 0 "SARIMA Model Prediction" # plot residual errors 'kde' The model summary and the predicted outcome are illustrated from the below plots. SARIMAX Model Summary Predicted vs Actual SARIMA In order to plot the trained, validation and predicted data under proper dates, we use the following code snippet. plt.figure(figsize=( , ))

index= pd.date_range(start=train_dates[ ], periods=len(train_dates), freq= )
valid_index = pd.date_range(start=valid_dates[ ], periods=len(valid_dates), freq= )

train_active =  pd.Series(train_ml[ ].values, index)
valid_active =  pd.Series(valid_ml[ ].values, valid_index)
pred_active =  pd.Series(prediction_ar, valid_index)

f, ax = plt.subplots( , , figsize=( , ))
plt.plot(train_active, marker= ,color= ,label = )
plt.plot(valid_active, marker= ,color= ,label = )
plt.plot(pred_active, marker= ,color= ,label = )

plt.legend()
plt.xlabel( )
plt.ylabel( )
plt.title( + stateName) 10 5 0 'D' 0 'D' 'Active Cases' 'Active Cases' 1 1 12 10 'o' 'blue' "Train Data Set" 'o' 'green' "Valid Data Set" 'o' 'red' "Predicted ARIMA/SARIMAX" "Date Time" 'Active Cases' "Active Cases ARIMA/SARIMAX Model Forecasting for state " Auto-Correlation Plots Differencing helps to make time-series stationary. In order to find a right order of differencing, we use the ACF plot and see that it reaches zero, fast and with a fairly static mean. With Covid-19 data for state Maharashtra, the ACF plot and the code snippet are given below. fig, axes = plt.subplots( , , sharex= , figsize=( , ))
axes[ , ].plot(df_per_State_features[ ]); axes[ , ].set_title( )
plot_acf(df_per_State_features[ ], ax=axes[ , ]) axes[ , ].plot(df_per_State_features[ ].diff()); axes[ , ].set_title( )
plot_acf(df_per_State_features[ ].diff().dropna(), ax=axes[ , ]) axes[ , ].plot(df_per_State_features[ ].diff().diff()); axes[ , ].set_title( )
plot_acf(df_per_State_features[ ].diff().diff().dropna(), ax=axes[ , ])

plt.show() # Original Series 3 2 True 10 5 0 0 'Active Cases' 0 0 'Original Series' 'Active Cases' 0 1 # 1st Differencing 1 0 'Active Cases' 1 0 '1st Order Differencing' 'Active Cases' 1 1 # 2nd Differencing 2 0 'Active Cases' 2 0 '2nd Order Differencing' 'Active Cases' 2 1 AUTO ARIMA MODELS This library automatically discovers the optimal order for an ARIMA model with stepwise execution of hyperparameters and parallel fitting of models. This objective of this library (auto_arima) is to identify the most optimal parameters for an ARIMA/SARIMA and return a fitted ARIMA model. It does not depend on the , but instead, it conducts differencing tests (i.e., Kwiatkowski–Phillips–Schmidt–Shin, Augmented Dickey-Fuller or Phillips–Perron) on its own to determine the order of differencing. PACF/Auto-Correlation (manual computation of differencing) One other difference from Stats Model’s Vanilla ARIMA is that in case the model does not converge, it raises ValueError to signify the re-evaluation of stationarity-inducing measures. For seasonality, SARIMA identifies the optimal P and Q hyperparameters after conducting the Canova-Hansen to determine the optimal order of seasonal differencing , so that the auto-regressive/moving average portions of the seasonal model are defined by start_P, start_Q, max_P, max_Q D The model is fitted within ranges of defined start_p order (or number of time lags for AR), max_p, start_q (order of the moving average MA), max_q ranges. Hence it can be formulated as : (p, q, d)x(P, Q, D) The best model is obtained by optimizing ( – Akaike Information Criterion, – Corrected Akaike Information Criterion, – Bayesian Information Criterion, – Hannan-Quinn Information Criterion, or “out of bag”–for validation scoring–respectively to a minimal value. information_criterion AIC AICC BIC HQIC It also operates with variables (just like state space methods/models) for predicting added features in the regression operation. exogenous Summary of AR with Auto-ARIMA The following code and figure depicts AR model with Auto ARIMA with start_p=0, start_q=2 (by default), max_p=5, max_q=0 model_ar= auto_arima(trainActiveCases,trace= , error_action= , start_p= ,start_q= ,max_p= ,max_q= ,
                   suppress_warnings= ,stepwise= ,seasonal= )
model_ar_fit = model_ar.fit(trainActiveCases)
prediction_ar=model_ar_fit.predict(len(validActiveCases))

print(model_ar_fit.summary())
plt.figure(figsize=( , )) residuals = pd.DataFrame(model_ar_fit.resid())
model_ar_fit.plot_diagnostics()
plt.show()
residuals.plot(kind= )
plt.show() True 'ignore' 0 0 5 0 True False False 10 10 # plot residual errors 'kde' Predicted vs Actual AR Auto-ARIMA Now let's look at the Moving Average model by changing (set to 5). Here's the code snippet for it. start_p (AutoRegressive Moving part - set to 0), start_q, max_p, max_q model_ma= auto_arima(trainActiveCases,trace= , error_action= , start_p= ,start_q= ,max_p= ,max_q= ,
                   suppress_warnings= ,stepwise= ,seasonal= )
model_ma_fit = model_ma.fit(trainActiveCases)
prediction_ma=model_ma_fit.predict(len(validActiveCases))

print(model_ma_fit.summary())
plt.figure(figsize=( , )) residuals = pd.DataFrame(model_ma_fit.resid())
model_ma_fit.plot_diagnostics()
plt.show()
residuals.plot(kind= )
plt.show() True 'ignore' 0 0 0 5 True False False 10 10 # plot residual errors 'kde' The hyper-parameter tuning for evaluating the lowest and best AIC with the residual components are given below. Summary of MA with Auto-ARIMA Predicted vs Actual MA Auto-ARIMA The Auto ARIMA is then used to model ARIMA by using the following parameters ------------------ as shown in the below code snippet. (start_p=1,start_q=1,max_p=5,max_q=5) model_arima= auto_arima(trainActiveCases,trace= , error_action= , start_p= ,start_q= ,max_p= ,max_q= ,
                   suppress_warnings= ,stepwise= ,seasonal= )
model_arima_fit = model_arima.fit(trainActiveCases)
prediction_arima=model_arima_fit.predict(len(validActiveCases))

print(model_arima_fit.summary())
plt.figure(figsize=( , )) residuals = pd.DataFrame(model_arima_fit.resid())
model_arima_fit.plot_diagnostics()
plt.show()
residuals.plot(kind= )
plt.show() True 'ignore' 1 1 5 5 True False False 10 10 # plot residual errors 'kde' Summary of ARIMA with Auto-ARIMA Predicted vs Actual Auto-ARIMA AUTO SARIMA MODEL The SARIMA model is tuned with to find the best model with the lowest AIC. The code snippet is given below, which does parameter tuning for Here and are tuned with values (1,0,1), (1,0,0),  (0,0,1)  and 12 respectively. hyper-parameters (p, q, d)x(P, Q, D) (P, Q, D) (m) The parameters are initialized as , . even , takes start_p=0,start_q=2,max_p=5,max_q=5,m=12 start_q set to 0 default value of 2. model_sarima= auto_arima(trainActiveCases,trace= , error_action= , 
                         start_p= ,start_q= ,max_p= ,max_q= ,m= ,
                   suppress_warnings= ,stepwise= ,seasonal= )
model_sarima_fit = model_sarima.fit(trainActiveCases)
prediction_sarima=model_sarima_fit.predict(len(validActiveCases))

print(model_sarima_fit.summary())
plt.figure(figsize=( , )) residuals = pd.DataFrame(model_sarima_fit.resid())
model_sarima_fit.plot_diagnostics()
plt.show()
residuals.plot(kind= )
plt.show() True 'ignore' 0 2 5 5 12 True True True 10 10 # plot residual errors 'kde' Summary of SARIMA with Auto-SARIMA Predicted vs Actual Auto-SARIMA Seasonal Differencing SARIMA We also try to explore seasonality with Covid-19 dataset by applying regular differencing, where the values are differenced from the previous season. In the equation of SARIMA(p,d,q)x(P,D,Q): P denotes Seasonal AR, D denote Seasonal Order of seasonal differencing and Q denotes Seasonal Moving Average, x is the frequency of the time series. If a model has defined seasonal patterns, we set D=1 for a given frequency ‘x’. For a SARIMA model, we keep the total differencing + D’ <=2 SELECTION OF BEST MODEL with AUTO-ARIMA The following figure illustrates how to select the best Auto ARIMA Model, compute its error metrics for any input time-series data. Comparative Study of RMSE Conclusion and Future Work More work needs to be done by studying the impact of the number of people quarantined, hospital beds, icu units, ventilators, rate of tests done by using Vector Autoregressive Moving-Average with eXogenous regressors (VARMAX). Evaluate the model performance by (multiplicative models) and the two different modules of StatsModel. sm.tsa.statespace.ExponentialSmoothing sm.tsa.ExponentialSmoothing (linear model), i.e. A detailed study of using different smoothing techniques with latency and accuracy to analyze the tradeoff between the two. Evaluate confidence intervals for forecasts, based on an assumption of Gaussian errors as provided by Analyzing the performance impact by concentrating initial values out of the objective function, particularly in situations when the seasonal periodicity is large. sm.tsa.statespace.ExponentialSmoothing With more number of COVID-19 cases evolving, the right order of differencing (if needed multiple orders of differencing) needs to be computed to arrive at a near-stationary series (ACF plot nearing zero). Here we have only computed RMSE, but an ideal situation all error metrics Mean Absolute Percentage Error (MAPE), Mean Error (ME), Mean Absolute Error (MAE), Mean Percentage Error (MPE), Root Mean Squared Error (RMSE), Lag 1 Autocorrelation of Error (ACF1), Correlation between the Actual and the Forecast (core), Min-Max Error (min-max) needs to be computed. The complete source code is available at https://github.com/sharmi1206/covid-19-analysis References https://www.statsmodels.org/dev/user-guide.html#time-series-analysis http://alkaline-ml.com/pmdarima/0.9.0/modules/generated/pyramid.arima.auto_arima.html https://www.machinelearningplus.com/time-series/arima-model-time-series-forecasting-python/ https://machinelearningmastery.com/how-to-grid-search-sarima-model-hyperparameters-for-time-series-forecasting-in-python/https://otexts.com/fpp2/arima-r.html

Initialized

Trace

Basic Understanding of ARIMA/SARIMA vs Auto ARIMA/SARIMA using Covid-19 Data Predictions

About Author

Comments

TOPICS

THIS ARTICLE WAS FEATURED IN

Related Stories

Anomaly Detection from Fetal ECG — A Case Study of IOT Anomaly Detection using GAN

10 Non-Influencers redefining Influence. Opportunity: $1 Million USD Funding to fight COVID-19

121 Stories To Learn About Coronavirus Impact On Business

278 Stories To Learn About Covid 19

196 Stories To Learn About Covid19

18 Hard Truths About Corona Recession

Anomaly Detection from Fetal ECG — A Case Study of IOT Anomaly Detection using GAN

10 Non-Influencers redefining Influence. Opportunity: $1 Million USD Funding to fight COVID-19

121 Stories To Learn About Coronavirus Impact On Business

278 Stories To Learn About Covid 19

196 Stories To Learn About Covid19

18 Hard Truths About Corona Recession

Light-Mode

Classic

Newspaper

Minty

Dark-Mode

Neon Noir

Minty

HN StartUps