<jupyter_start><jupyter_text>ForecastThe stationarity of the time series for women employmenmt data was investigated in Notebook "Timeseries Analysis Stationarity".In this notebook, some forecast models are used for prediction.<jupyter_code>import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import statsmodels.api as sm
from statsmodels.tsa.holtwinters import SimpleExpSmoothing, Holt, ExponentialSmoothing
from scipy import signal
from sklearn.metrics import mean_absolute_error, mean_squared_error
import calendar

%matplotlib inline<jupyter_output><empty_output><jupyter_text>Utility functions<jupyter_code>def TrainTestSplit(Df, TrainPcnt = 0.75, FcstH = 24):
    """
    Split data into training and test sets
    Df: pandas dataframe in the form of [index, values]
    TrainPcnt: percentage of data allotted for training
    FcstH: Forcast horizon
    """
    # Train-Test Split

    Len = len(Df)
    TrainLen = int(TrainPcnt * Len)
    TrainDf = Df[:TrainLen]
    TestDf = Df[TrainLen:TrainLen + FcstH]
    return TrainDf, TestDf    

def PlotDiagnostics(Residu):
    """
    Plot diagnostics.
    """
    fig = plt.figure(figsize=(12,8))
    grid = plt.GridSpec(2,2, wspace=0.4, hspace=0.4)
    
    ax_resid = plt.subplot(grid[0,0])
    ax_hist = plt.subplot(grid[0,1])
    ax_qq = plt.subplot(grid[1,0])
    ax_acf = plt.subplot(grid[1,1])
    
    ax_resid.plot(Residu)
    ax_resid.set_xlabel('index')
    ax_resid.set_ylabel('residual')
    
    sns.distplot(Residu, hist=True, kde=True, 
             color = 'darkblue', 
             hist_kws={'edgecolor':'black'},
             kde_kws={'linewidth': 4}, ax=ax_hist)
    ax_hist.set_xlabel('residu')
    ax_hist.set_ylabel('density')
    
    sm.qqplot(Residu.values.flatten(), fit=True, line='45', ax=ax_qq)
    
    sm.tsa.graphics.plot_acf(Residu, ax=ax_acf)
def PlotForecast(TrainDf, TestDf, FittedDf, ForecastDf, ModelName=''):
    """
    Plot data, fitted values and forecast values
    MdlFit: Fitted Model e.g. MdlFit = Mdl.fit()
    """
    plt.figure(figsize=(12,16))
    grid = plt.GridSpec(3,1, wspace=0.4, hspace=0.4)
    ax_data = plt.subplot(grid[0,0])
    ax_fcst = plt.subplot(grid[1,0])
    
    ax_data.plot(TrainDf,'b',label='employment data')
    ax_data.plot(FittedDf,'r', label='Fitted value '+ModelName)
    ax_data.set_ylabel('vlaue')
    ax_data.set_title('Fitted exponential smoothing, Mode: '+ModelName)
    ax_data.legend()
    
    LenTrainToPlot=len(TrainDf)//10
    ax_fcst.plot(TrainDf[-LenTrainToPlot:],'b', label='Train data')
    ax_fcst.plot(FittedDf[-LenTrainToPlot:],'r', label='fitted')
    ax_fcst.plot(TestDf,'g', label='Test data')
    ax_fcst.plot(ForecastDf,'k--', label='Forecast '+ModelName)
    
    ax_fcst.set_ylabel('vlaue')
    ax_fcst.set_title('Forecast exponential smoothing, Mode: '+ModelName)
    ax_fcst.legend()
def Performance(TrainDf, TestDf, FittedDf, ForecastDf, ModelName=''):
    
def PlotSeasonal(TS, Period='year'):
    """
    Plot seasonal data of one time series
    TS: pandas series
    """
    # Colors
    Colors = ['tab:red', 'tab:blue', 'tab:green', 'tab:orange', 'tab:brown', 'tab:grey', 'tab:pink', 'tab:olive', 'deeppink', 'steelblue', 'firebrick', 'mediumseagreen']      
    NumColors = len(Colors)
    plt.figure(figsize=(16,10), dpi= 80)
    ax = plt.subplot()
    # group by year
    TS.sort_index(inplace=True)
    gb = TS.groupby(TS.index.year)
    NumGps = gb.ngroups
    
    # plot
    Months = np.array(calendar.month_abbr)
    for i,yr in enumerate(gb.groups.keys()):
        ts = gb.get_group(yr)
        ts.index = ts.index.month
        ax.plot(ts,Colors[i%NumColors], label=yr)
        ax.text(ts.index[-1]+0.1, ts.values[-1], yr, fontsize=12, color = Colors[i%NumColors] )
    ax.set_xticks(range(13))
    ax.set_xticklabels(Months)
    plt.xlabel('Months')
    plt.ylabel('Series Value')

def NaivePeakTroughFinder(TS, bPlot=False):
    """
    Find peaks and troughs using differencing
    """
    TsVal = TS.values
    Diff1 = np.sign(np.diff(TsVal))
    Diff2 = np.diff(Diff1)
    PkIntIndx = np.where(Diff2==-2)[0] + 1
    TrIntIndx = np.where(Diff2== 2)[0] + 1
    
    PkIndx = TS.index[PkIntIndx]
    TrIndx = TS.index[TrIntIndx]
    
    # plot
    if(bPlot):
        plt.figure(figsize=(14,4))
        ax=plt.subplot()
        ax.plot(TS,label='time series')
        ax.plot(PkIndx, TS[PkIndx], '^r', label='Peaks')
        ax.plot(TrIndx, TS[TrIndx], 'vg', label='Trough')
    
    return PkIndx, TrIndx<jupyter_output><empty_output><jupyter_text>Women employment rate in the USdata source: https://fred.stlouisfed.org/<jupyter_code>WmnEmpRate = pd.read_csv('women_emp_rate_25_52.csv') 
WmnEmpRate.head()
WmnEmpRate.info()
WmnEmpRate.columns=['date', 'rate']  # rename columns
WmnEmpRate['date'] = pd.to_datetime(WmnEmpRate['date']) # datetime column type
WmnEmpRate.set_index('date', inplace=True)
WmnEmpRate.head()<jupyter_output><empty_output><jupyter_text>take a look at the data and the rolling mean with window = 12 months<jupyter_code>plt.figure(figsize=(16,4))
plt.plot(WmnEmpRate,label='emp. rate')
plt.plot(WmnEmpRate.rolling(12).mean(), 'r',label='rolling mean 12 mo' )
plt.xlabel('date')
plt.ylabel('percent')
plt.title('women employment rate')
plt.legend()<jupyter_output><empty_output><jupyter_text>There is clearly a trand in the mean, but the variance seems to be stable. There are many ways to remove the trend in the data.For time series for which the waveform morphology is important, it's better to use a technique with minimal distorsion to the waveform; for example, a zero distortion high pass FIR filtering. Other methods such as differencing are as effective for manyapplications.  Forecasting: Exponential Smoothing (ES)Exponential smoothing are simple and fast forecasting algorithms. Table below list the caplability of theseforecasting schemes in handeling time series with various characteristics.Algorithm| stationary | trend |seasonality|--|--|--|--|Simple ES| y | - |-|Holt's ES| y | y |-|Winter's ES| y | y |y|In simple ES,$$\hat{y}_{t+h|t} = L_t$$$$L_t = \alpha Y_{t-1} + (1-\alpha) L_{t-1}$$the level $L_t$ is updated as a weighted average of the value at $t$ and the previous levelwhere $0 \le \alpha \le 1$.Holt improved the forecast for time series with trend $$\hat{y}_{t+h|t} = L_t + h\, T_t$$$$L_t = \alpha Y_t + (1-\alpha) (L_{t-1} + T_{t-1})$$$$T_t = \beta (L_t - L_{t-1}) + (1-\beta) T_{t-1}$$where trend $T_t$ is updated as a weigthed avergae of the current and previous With Winter's improvement the seasonality can be incorporated into model.  Train-Test Split<jupyter_code>TrainPcnt = 0.75
FcstH = 24
TrainDf, TestDf = TrainTestSplit(WmnEmpRate, TrainPcnt, FcstH)
print("Train len: {0:d}, Test len: {1:d}".format(len(TrainDf), len(TestDf)))<jupyter_output>Train len: 540, Test len: 24<jupyter_text>Simple Exponential Smoothing (SES)<jupyter_code>SESMdl = SimpleExpSmoothing(TrainDf)
SESFit = SESMdl.fit(smoothing_level=.2)
#PredSES = SESFit.forecast(FcstH)
PredSES = SESFit.predict(start=TestDf.index[0],end=TestDf.index[-1])
SESFit.summary()
PlotForecast(TrainDf, TestDf, SESFit.fittedvalues, PredSES, 'SES')
PlotDiagnostics(SESFit.resid)<jupyter_output><empty_output><jupyter_text>Holt's Exponentail Smoothing HES<jupyter_code>HESMdl = Holt(TrainDf)
HESFit = HESMdl.fit(smoothing_level=0.2, smoothing_slope=0.2)
PredHES = HESFit.predict(start=TestDf.index[0],end=TestDf.index[-1])
PlotForecast(TrainDf, TestDf, HESFit.fittedvalues, PredHES, 'HES')
PlotDiagnostics(HESFit.resid)
HESFit.summary()<jupyter_output><empty_output><jupyter_text>Winter's Exponential Smoothing WES<jupyter_code>WESMdl = ExponentialSmoothing(TrainDf, trend='add', seasonal='add', seasonal_periods=12, damped=False)
WESFit = WESMdl.fit(smoothing_level=0.2, smoothing_slope=0.2, smoothing_seasonal=0.3)
PredWES = WESFit.predict(start=TestDf.index[0],end=TestDf.index[-1])
PlotForecast(TrainDf, TestDf, WESFit.fittedvalues, PredWES, 'HES')
PlotDiagnostics(WESFit.resid)
WESFit.summary()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13805090006207324, 'fraction_numerical': 0.015394165114835505, 'mean_word_length': 4.602225312934632, 'pattern_counts': {'":': 0, '<': 32, '<?xml version=': 0, '>': 32, 'https://': 1, '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': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '67297', 'license': 'MIT', 'n_tokens_mistral': 3182, 'n_tokens_neox': 2999, 'n_words': 627, 'path': 'Forecast_Exponentail_Smoothing.ipynb', 'repo_name': 'matti-logiko/Time-Series-Forecast-Women-Employment-Data', 'size': '8055'}