<jupyter_start><jupyter_text>Handling time series data in Pandas few operations<jupyter_code>pd.datetime.strptime('1949-01', '%Y-%m')
from dateutil.parser import parse
parse('1949-01')
data.index
data.loc['1949-01-23']
data['Passengers']
data.loc['1949-01-10':'1949-05-10']<jupyter_output><empty_output><jupyter_text>Check Stationarity of a Time SeriesA TS is said to be stationary if its statistical properties such as mean, variance remain constant over time. But why is it important? Most of the TS models work on the assumption that the TS is stationary. Intuitively, we can state that if a TS has a particular behaviour over time, there is a very high probability that it will follow the same in the future. Also, the theories related to stationary series are more mature and easier to implement as compared to non-stationary series.Stationarity is defined using very strict criterion. However, for practical purposes we can assume the series to be stationary if it has constant statistical properties over time, ie. the following:* constant mean* constant variance* an autocovariance that does not depend on time<jupyter_code>data.plot() # plotting with pandas
plt.plot(data)<jupyter_output><empty_output><jupyter_text>Plotting Rolling Statistics:We can plot the moving average or moving variance and see if it varies with time. By moving average/variance I mean that at any instant ‘t’, we’ll take the average/variance of the last year, i.e. last 12 months. But again this is more of a visual technique<jupyter_code>data.Passengers.plot()
data.Passengers.rolling(12).mean().plot(label='rolling mean')
data.Passengers.rolling(12).std().plot(label='rolling std')
plt.legend()
plt.title('Rolling Mean & Standard Deviation')
data.Passengers.rolling(12).mean()<jupyter_output><empty_output><jupyter_text>Perform Dickey-Fuller test:Dickey-Fuller Test: This is one of the statistical tests for checking stationarity. Here the null hypothesis is that the TS is non-stationary. The test results comprise of a Test Statistic and some Critical Values for difference confidence levels. If the ‘Test Statistic’ is less than the ‘Critical Value’, we can reject the null hypothesis and say that the series is stationary. Refer this article for details.<jupyter_code>from statsmodels.tsa.stattools import adfuller
data_series=(data.values).ravel() # in full
print('Results of Dickey-Fuller Test:')
dftest = adfuller(data_series, autolag='AIC')
dfoutput = pd.Series(dftest[0:4], index=['Test Statistic','p-value','#Lags Used','Number of Observations Used'])
for key,value in dftest[4].items():
    dfoutput['Critical Value (%s)'%key] = value
print(dfoutput)<jupyter_output>Results of Dickey-Fuller Test:
Test Statistic                   0.815369
p-value                          0.991880
#Lags Used                      13.000000
Number of Observations Used    130.000000
Critical Value (1%)             -3.481682
Critical Value (5%)             -2.884042
Critical Value (10%)            -2.578770
dtype: float64<jupyter_text>Though the variation in standard deviation is small, mean is clearly increasing with time and this is not a stationary series. Also, the test statistic is way more than the critical values. Note that the signed values should be compared and not the absolute values.Next, we’ll discuss the techniques that can be used to take this TS towards stationarity.  How to make a Time Series Stationary?Though stationarity assumption is taken in many TS models, almost none of practical time series are stationary. So statisticians have figured out ways to make series stationary, which we’ll discuss now. Actually, its almost impossible to make a series perfectly stationary, but we try to take it as close as possible.Lets understand what is making a TS non-stationary. There are 2 major reasons behind non-stationarity of a TS:*  Trend – varying mean over time. For eg, in this case we saw that on average, the number of passengers was growing over time.*  Seasonality – variations at specific time-frames. eg people might have a tendency to buy cars in a particular month because of pay increment or festivals. STEPS* Step:1> The underlying principle is to model or estimate the trend and seasonality in the data and remove those from the data to get a stationary series. * Step:2> Statistical forecasting techniques can be implemented on this series. The final step would be to convert the forecasted values into the original scale by applying trend and seasonality constraints back.Note: Here we will be discussing a number of methods. Some might work well in this case and others might not. But the idea is to get a hang of all the methods and not focus on just the problem at hand.Working on the trend part:  Estimating & Eliminating TrendOne of the first tricks to reduce trend can be transformation. For example, in this case we can clearly see that the there is a significant positive trend. So we can apply transformation which penalize higher values more than smaller values. These can be taking a log, square root, cube root, etc. Lets take a log transform here for simplicity:<jupyter_code>## Log transformation
data_log = np.log(data)
data_log.plot()<jupyter_output><empty_output><jupyter_text>It is easy to see a forward trend in the data. But its not very intuitive in presence of noise. So we can use some techniques to estimate or model this trend and then remove it from the data. Most commonly used techniques are:  1. Aggregation – taking average for a time period like monthly/weekly averages  2. Smoothing – taking rolling averages  3. Polynomial Fitting – fit a regression modelWe will use smoothing here and you should try other techniques as well which might work out for other problems. Smoothing refers to taking rolling estimates, i.e. considering the past few instances. There are can be various ways,  two are discussed here.  Moving average<jupyter_code>data_log.Passengers.plot()
data_log.Passengers.rolling(12).mean().plot(label='rolling mean')
plt.legend()
plt.title('Rolling Mean after log tranformation')<jupyter_output><empty_output><jupyter_text>Subtracting  the rolling mean  from the original series. Note that since we are taking average of last 12 values, rolling mean is not defined for first 11 values.<jupyter_code>ra=data_log.Passengers.rolling(12).mean()
data_diff = data_log.Passengers- ra
data_diff.dropna(inplace=True)
data_diff
data_diff.plot()
data_diff.rolling(12).mean().plot(label='rolling mean')
data_diff.rolling(12).std().plot(label='rolling std')
plt.legend()
plt.title('Rolling Mean & Standard Deviation')
### Again checking Dickey-Fuller Test
print('Results of Dickey-Fuller Test:')
dftest = adfuller(data_diff, autolag='AIC')
dfoutput = pd.Series(dftest[0:4], index=['Test Statistic','p-value','#Lags Used','Number of Observations Used'])
for key,value in dftest[4].items():
    dfoutput['Critical Value (%s)'%key] = value
print(dfoutput)<jupyter_output>Results of Dickey-Fuller Test:
Test Statistic                  -3.162908
p-value                          0.022235
#Lags Used                      13.000000
Number of Observations Used    119.000000
Critical Value (1%)             -3.486535
Critical Value (5%)             -2.886151
Critical Value (10%)            -2.579896
dtype: float64<jupyter_text>This looks like a much better series. The rolling values appear to be varying slightly but there is no specific trend. Also, the test statistic is smaller than the 5% critical values so we can say with 95% confidence that this is a stationary series.  Exponentially Weighted Moving Average  In above approach time-period has to be strictly defined. In this case we can take yearly averages but in complex situations like forecasting a stock price, its difficult to come up with a number. So we take a ‘weighted moving average’ where more recent values are given a higher weight. There can be many technique for assigning weights. A popular one is exponentially weighted moving average where weights are assigned to all the previous values with a decay factor.This can be implemented in Pandas as:<jupyter_code>ewa = data_log.ewm(span=12).mean()
plt.plot(data_log,label='Log transformed ori. data')
plt.plot(data_log.rolling(12).mean(),label='rolling mean')
plt.plot(ewa,label='exp. wtd moving average')
plt.legend()
data_diff1 = data_log- ewa
data_diff1
plt.plot(data_diff1)
plt.plot(data_diff1.rolling(12).mean(),label='rolling mean')
plt.plot(data_diff1.rolling(12).std(),label='rolling std')
plt.legend()
plt.title('Rolling Mean & Standard Deviation')
data_diff1.Passengers ## convverting data into series
## Again checking Dickey-Fuller test
print('Results of Dickey-Fuller Test:')
dftest = adfuller(data_diff1.Passengers, autolag='AIC')
dfoutput = pd.Series(dftest[0:4], index=['Test Statistic','p-value','#Lags Used','Number of Observations Used'])
for key,value in dftest[4].items():
    dfoutput['Critical Value (%s)'%key] = value
print(dfoutput)<jupyter_output>Results of Dickey-Fuller Test:
Test Statistic                  -3.310542
p-value                          0.014414
#Lags Used                      13.000000
Number of Observations Used    130.000000
Critical Value (1%)             -3.481682
Critical Value (5%)             -2.884042
Critical Value (10%)            -2.578770
dtype: float64<jupyter_text>Eliminating Trend and SeasonalityThe simple trend reduction techniques discussed before don’t work in all cases, particularly the ones with high seasonality. Lets discuss two ways of removing trend and seasonality:* Differencing – taking the differece with a particular time lag* Decomposition – modeling both trend and seasonality and removing them from the model. Differencing :>In this technique, we take the difference of the observation at a particular instant with that at the previous instant. This mostly works well in improving stationarity. First order differencing :difference(t) = observation(t) - observation(t-1), can be done in Pandas as given below ><jupyter_code># create a differenced series # input data must be in seris 
def difference(dataset, interval=1):
    diff = list()
    for i in range(interval, len(dataset)):
        value = dataset[i] - dataset[i - interval]
        diff.append(value)
    return pd.Series(diff)
plt.plot(difference(data_log.values))
data_log_diff = data_log.diff() # Default period is taken 1
plt.plot(data_log_diff)
plt.show()
data_log_diff.dropna(inplace=True)
ts_log_diff = data_log - data_log.shift(1)
plt.plot(ts_log_diff)
## Again plot rolling mean and rolling std and chec Dickey fuller test<jupyter_output><empty_output><jupyter_text>Breaking a time series: * Systematic: Components of the time series that have consistency or recurrence and can be described and modeled are viz.evel, trend, seasonality.* Non-Systematic: Components of the time series that cannot be directly modeled is noiseThese components are defined as follows:* Level: The average value in the series.* Trend: The increasing or decreasing value in the series.* Seasonality: The repeating short-term cycle in the series.* Noise: The random variation in the series.>So a time series is thought to be an aggregate or combination of these four components.>All series have a level and noise. The trend and seasonality components are optional.>It is helpful to think of the components as combining either additively or multiplicatively as given by relation below:* y(t) = Level + Trend + Seasonality + Noise* y(t) = Level * Trend * Seasonality * Noise<jupyter_code>#Multiplicative decomposition
from statsmodels.tsa.seasonal import seasonal_decompose
result = seasonal_decompose(data_log, model='multiplicative',freq=12)
result.plot()
plt.show()
trend = result.trend
seasonal = result.seasonal
residual = result.resid
residual.dropna(inplace=True)
residual<jupyter_output><empty_output><jupyter_text>Again check rolling mean and std and Dickey fuller test:You will find that test statistic is significantly lower than the 1% critical value. So this TS is very close to stationary. All above are the techniques for making time series data stationary* A strictly stationary series with no dependence among the values. This is the easy case wherein we can model the residuals as noise. But this is very rare.* A series with significant dependence among values. In this case we need to use some statistical models like ARIMA to forecast the data. ARIMA>ARIMA: Auto-Regressive Integrated Moving Averages>The ARIMA forecasting for a stationary time series is nothing but a linear (like a linear regression) equation. The predictors depend on the parameters (p,d,q) of the ARIMA model:![Arima%20explanation.PNG](attachment:Arima%20explanation.PNG)  WE will work on differenced data :> data_log_diffDifferencing is a popular technique and relatively easier to add noise and seasonality back into predicted residuals in this case.<jupyter_code>#ACF and PACF plots:
from statsmodels.tsa.stattools import acf, pacf
lag_acf = acf(data_log_diff, nlags=12)
lag_pacf = pacf(data_log_diff, nlags=12)
#Plot ACF: 
plt.subplot(121) 
plt.plot(lag_acf)
plt.axhline(y=0,linestyle='--',color='gray')
plt.axhline(y=-1.96/np.sqrt(len(data_log_diff)),linestyle='--',color='gray')
plt.axhline(y=1.96/np.sqrt(len(data_log_diff)),linestyle='--',color='gray')
plt.title('Autocorrelation Function')
#Plot PACF:
plt.subplot(122)
plt.plot(lag_pacf)
plt.axhline(y=0,linestyle='--',color='gray')
plt.axhline(y=-1.96/np.sqrt(len(data_log_diff)),linestyle='--',color='gray')
plt.axhline(y=1.96/np.sqrt(len(data_log_diff)),linestyle='--',color='gray')
plt.title('Partial Autocorrelation Function')
plt.tight_layout()<jupyter_output><empty_output><jupyter_text>In this plot, the two dotted lines on either sides of 0 are the confidence interevals. These can be used to determine the ‘p’ and ‘q’ values as:![acf%20pcf.PNG](attachment:acf%20pcf.PNG)<jupyter_code>from statsmodels.graphics.tsaplots import plot_acf
plot_acf(data_log_diff, lags=12)
plt.show()
from statsmodels.graphics.tsaplots import plot_pacf
plot_pacf(data_log_diff, lags=12,use_vlines=True, marker='o')
plt.show()
from statsmodels.tsa.arima_model import ARIMA
#The p,d,q values can be specified using the 
#order argument of ARIMA which take a tuple (p,d,q).
#Let model the 3 cases:
data_log_diff.Passengers
## AR Model
model = ARIMA(data_log, order=(2, 1, 0))  
results_AR = model.fit(disp=-1)
predict1=results_AR.fittedvalues
plt.plot(data_log_diff)
plt.plot(predict1, color='red')
plt.title('RSS: %.4f'% sum((predict1-(data_log_diff.Passengers))**2))
## MA Model
model = ARIMA(data_log, order=(0, 1, 2))  
results_MA = model.fit(disp=-1)
predict2=results_MA.fittedvalues
plt.plot(data_log_diff)
plt.plot(predict2, color='red')
plt.title('RSS: %.4f'% sum((predict2-(data_log_diff.Passengers))**2))
## Combined Model
model = ARIMA(data_log, order=(2, 1, 2))  
results_ARIMA = model.fit(disp=-1)  
predict3=results_ARIMA.fittedvalues ## one oder differnced values
### Thus comparing with differenced data 'data_log_diff'
plt.plot(data_log_diff)
plt.plot(predict3, color='red')
plt.title('RSS: %.4f'% sum((predict3-(data_log_diff.Passengers))**2))<jupyter_output>C:\Users\hP\Anaconda3\lib\site-packages\statsmodels\tsa\base\tsa_model.py:218: ValueWarning: A date index has been provided, but it has no associated frequency information and so will be ignored when e.g. forecasting.
  ' ignored when e.g. forecasting.', ValueWarning)
C:\Users\hP\Anaconda3\lib\site-packages\statsmodels\tsa\base\tsa_model.py:218: ValueWarning: A date index has been provided, but it has no associated frequency information and so will be ignored when e.g. forecasting.
  ' ignored when e.g. forecasting.', ValueWarning)<jupyter_text>Converting Back the data into original scale:<jupyter_code>predictions_ARIMA_diff = predict3.copy()
predictions_ARIMA_diff
predictions_ARIMA_diff_cumsum = predictions_ARIMA_diff.cumsum()
predictions_ARIMA_diff_cumsum<jupyter_output><empty_output><jupyter_text>Adding first value of modeled data (i.e. data_log)  to each value of 'predicitons_ARIMA_diff_cumsum'<jupyter_code>data_log
a=np.ones(len(data_log))*data_log.values[0][0]
p_first_value =pd.Series(a,index=data_log.index)
# finally
predictions_ARIMA_log=p_first_value.add(predictions_ARIMA_diff_cumsum,fill_value=0)
predictions_ARIMA_log
predictions_ARIMA = np.exp(predictions_ARIMA_log)
plt.plot(data,label='Initial')
plt.plot(predictions_ARIMA, label='Predicted')
plt.title('RMSE: %.4f'% np.sqrt(sum((predictions_ARIMA-data.Passengers)**2)/len(data)))
plt.legend()<jupyter_output><empty_output><jupyter_text>Forcasting for future one time step ahead<jupyter_code># data_log has index but it is monthly data we should convert it to DatetimeIndex and pass a frequency by to_period:
data_log.index = pd.DatetimeIndex(data_log.index).to_period('M')
model1 = ARIMA(data_log, order=(2, 1, 2))  
model_fit1 = model.fit()  
output1=model_fit1.forecast() ## one time step ahead by default, we can pass number of time step in future we want to predict 
output1[0]
## Note forecast gives data in same scale as that of input data and we dont have to
## integrate it by applying cumsum and addition of first value
one_step_head=np.exp(output1[0]) ## conveting back to its original scale by taking exponential
one_step_head
# The given last value in the data set
data.iloc[-1]<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.0878327654905584, 'fraction_numerical': 0.02084656695732517, 'mean_word_length': 5.211621233859398, 'pattern_counts': {'":': 0, '<': 57, '<?xml version=': 0, '>': 67, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '15', 'dir': 'jupyter-structured-clean-dedup', 'id': '308421', 'license': 'MIT', 'n_tokens_mistral': 5420, 'n_tokens_neox': 5043, 'n_words': 1936, 'path': 'Air Passengers Forecasting.ipynb', 'repo_name': 'Rami-RK/Time_Series', 'size': '17317'}