<jupyter_start><jupyter_text>"Understanding the timeseries aspect on the air quality of Delhi"> "Objective is to understand time series modeling and in turn forecast air quality of Delhi"- toc:true- branch: master- badges: true- comments: true- categories: [pollution, time-series, pandas, visualization, prophet] It is stated very normally that the air quality is pretty worse during/after diwali period. In this data exercise I would be experimenting with the time series aspect of the air quality using facebook's api "Prophet". I would like to understand how weekends, festival season and general seasonality impact the air quality index for Delhi city<jupyter_code>#collapse_show
#importing liberaries
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import os
import warnings
from fbprophet import Prophet 
import statsmodels.api as sm
warnings.filterwarnings("ignore")
%matplotlib inline
#collapse_show
#reading the datase
df=pd.read_csv('city_day.csv')
#extracting the data for date and AQI
dfd=df[df['City']=='Delhi'][['Date','AQI']].reset_index()
dfd=dfd.drop(labels='index', axis=1)
#collapse_show
#inspecting the extracted dataset
#this is the entire city data for Delhi
dfd.info()<jupyter_output><class 'pandas.core.frame.DataFrame'>
RangeIndex: 2009 entries, 0 to 2008
Data columns (total 2 columns):
 #   Column  Non-Null Count  Dtype  
---  ------  --------------  -----  
 0   Date    2009 non-null   object 
 1   AQI     1999 non-null   float64
dtypes: float64(1), object(1)
memory usage: 31.5+ KB<jupyter_text>This shows that the data is comprised of two coloumns one Date and the other is AQI. We need to convert the Date which is an object into Datetime object.<jupyter_code>#collapse_show
#converting the same to data object
dfd['Date']=pd.to_datetime(df['Date'])
#collapse_show
dfd.info() #As can be seen here that the date time is converted from object to Datetime
#collapse_show
#plotting the time series data using matplotlib 
dfd.plot(x='Date', y='AQI',kind='line',
        title="Primary Graph-AQI vs year", grid= True, figsize=(20, 5))
#collapse_show
#Checking for missing datas
dfd.isna().sum()<jupyter_output><empty_output><jupyter_text>As we can see we have 10 missing values in the dataset<jupyter_code>#collapse_show
#imputing the missing values
#Checking which fields have missingg values
dfd[dfd['AQI'].isna()]<jupyter_output><empty_output><jupyter_text>Here, we can see that majorily the data is abscent from 22nd to 29th Aug for 2017. Correcting the same.<jupyter_code>#collapse_show
dfd['AQI']=dfd['AQI'].fillna(dfd['AQI'].rolling(4, min_periods=1).mean())
#collapse_show
dfd['AQI'].isna().sum()<jupyter_output><empty_output><jupyter_text>Exploratory Data AnalysisWe will first see the trend of the data with respect to different date range.<jupyter_code>#collapse_show
#creating a function to extract different date parameters
def date_feature(df, label=None):
    df= df.copy()
    
    df['date']= df['Date']
    df['month']=df['Date'].dt.strftime('%B')
    df['year']=df['Date'].dt.strftime('%Y')
    df['dayofweek']=df['Date'].dt.strftime('%A')
    df['quarter']=df['Date'].dt.quarter
    df['dayofyear']=df['Date'].dt.dayofyear
    df['dayofmonth']=df['Date'].dt.day
    df['weekofyear']=df['Date'].dt.weekofyear
    X= df[['date','month','year','dayofweek','quarter','dayofyear','dayofmonth','weekofyear']]
    if label:
        y= df[label]
        return X, y
    return X
X, y= date_feature(dfd, label='AQI')
#collapse_show
df_new=pd.concat([X,y], axis=1)
df_new.head()
#collapse_show
# Plotting different features to  see differnet trends
plt.figure(figsize=(12,5))
pallette= sns.color_palette("mako_r", 4)
a=sns.barplot(data=df_new, x='month', y='AQI', hue='year')
a.set_title("Month wise trend for AQI")
plt.legend(loc='best')
plt.show()<jupyter_output><empty_output><jupyter_text>We can see from this graph that the AQI is best in the summser season and goes to worst in winter season (esp bad in the month of Nov. due to Parali (Stubble buring and Diwali)<jupyter_code>#collapse_show
plt.figure(figsize=(12,5))
pallette= sns.color_palette("mako_r", 4)
df_agg=df_new.groupby('year')['AQI'].mean().reset_index()
a=sns.barplot(data=df_agg, x='year', y='AQI', hue='year')
a.set_title("year wise trend for AQI")
plt.legend(loc='best')
plt.show()<jupyter_output><empty_output><jupyter_text>In general the air quality is getting better for Delhi.<jupyter_code>#collapse_show
plt.figure(figsize=(12,5))
pallette= sns.color_palette("mako_r", 4)
df_agg=df_new.groupby('month')['AQI'].mean().sort_values(ascending=False).reset_index()
a=sns.barplot(data=df_agg, x='month', y='AQI')
a.set_title("month wise trend for AQI")
#plt.legend(loc='best')
plt.show()<jupyter_output><empty_output><jupyter_text>As it can be seen from the graph above that the worst air quality is in the month of November and the best air quality is in the month of July, August. This could be also due to the fact that  genrally in these months monsoon arrives in Delhi resulting in better AQI vis a vis remaining months.<jupyter_code>#collapse_show
plt.figure(figsize=(12,5))
pallette= sns.color_palette("mako_r", 4)
df_agg=df_new.groupby('dayofweek')['AQI'].mean().sort_values(ascending=False).reset_index()
a=sns.barplot(data=df_agg, x='dayofweek', y='AQI')
a.set_title("Day of week wise trend for AQI")
#plt.legend(loc='best')
plt.show()<jupyter_output><empty_output><jupyter_text>Not much of a difference here. AQI is marginally better on Monday/Tuesday and Worst on Friday.<jupyter_code>#collapse_show
plt.figure(figsize=(12,5))
pallette= sns.color_palette("mako_r", 4)
df_agg=df_new.groupby('quarter')['AQI'].median().sort_values(ascending=False).reset_index()
a=sns.barplot(data=df_agg, x='quarter', y='AQI')
a.set_title("month wise trend for AQI")
#plt.legend(loc='best')
plt.show()<jupyter_output><empty_output><jupyter_text>This further reiterates the point that the fourth quarter has most severe values for AQI<jupyter_code>#collapse_show
plt.figure(figsize=(12,5))
pallette= sns.color_palette("mako_r", 4)
df_agg=df_new.groupby('dayofmonth')['AQI'].median().sort_values(ascending=False).reset_index()
a=sns.barplot(data=df_agg, x='dayofmonth', y='AQI')
a.set_title("days of month wise trend for AQI")
#plt.legend(loc='best')
plt.show()<jupyter_output><empty_output><jupyter_text>No perticular trend visible here.<jupyter_code>#collapse_show
plt.figure(figsize=(12,5))
pallette= sns.color_palette("mako_r", 4)
df_agg=df_new.groupby('weekofyear')['AQI'].median().sort_values(ascending=False).reset_index()
a=sns.barplot(data=df_agg, x='weekofyear', y='AQI')
a.set_title("week of the year wise trend for AQI")
#plt.legend(loc='best')
plt.show()<jupyter_output><empty_output><jupyter_text>The trend follows closely with the month trend ----  Modeling the data usin fb prophet We would be using the dfd dataframe. Also, dividing the data into training and test set.<jupyter_code>#collapse_show
f"the total number of days considered are:{dfd['Date'].count()}"<jupyter_output><empty_output><jupyter_text>As of now we are doing univariate analysis. Also, prophet expects the date to be specified as "ds" and output to be specified as "y"<jupyter_code>#collapse_show
#remaining the column name

dfd.columns=['ds','y']
dfd.columns
#collapse_show
#using 70% of the data as test data

loc1=int(0.7*2009)
dfd.loc[loc1]
mask1=dfd['ds']<='2018-11-07'
mask2=dfd['ds']>'2018-11-07'
X_train=dfd.loc[mask1]
X_test=dfd.loc[mask2]

f"the training size is: {X_train.shape}",f"the test size is: {X_test.shape}"
#collapse_show
#plotiing the train and test data in a single plot
fig,ax= plt.subplots(figsize=(12,5))
X_train.plot(kind='line', x='ds', y='y', color='green',ax=ax)
X_test.plot(kind='line', x='ds', y='y', color='red',ax=ax)
plt.title("Train vs Test distribution")
plt.grid(True)<jupyter_output><empty_output><jupyter_text>Around 70% of the data is being used for training and remaining 30% is being used as test data  Now we need a function to determine how far the prediction is from the actual value<jupyter_code>#collapse_show
#creating the error function

def mean_absolute_per_error(y_true,y_pred):
    y_true, y_pred= np.array(y_true), np.array(y_pred)
    return np.mean(np.abs(((y_true-y_pred)/y_true)*100))<jupyter_output><empty_output><jupyter_text>Creating a simple prophet model<jupyter_code>#collapse_show
model=Prophet()
model.fit(X_train)
#collapse_show
#doing a simple prediction
future=model.make_future_dataframe(periods=602, freq='D')
forecast=model.predict(future)
#collapse_show
forecast[['ds','yhat','yhat_lower','yhat_upper']].tail()<jupyter_output><empty_output><jupyter_text>The values as predicted are mentioned here. now comparing these values with the actual values<jupyter_code>#collapse_show
X_train.tail()<jupyter_output><empty_output><jupyter_text>As we can see there is a decent difference in prediction numbers.<jupyter_code>#collapse_show
#plotting the components of the model

fig=model.plot_components(forecast)<jupyter_output><empty_output><jupyter_text>Comparing actual Data with  forecast<jupyter_code>#collapse_show
X_test_forecast=model.predict(X_test)
X_test_forecast.tail()
#collapse_show
#comparing the test AQI with forecast AQI

f,ax= plt.subplots(figsize=(12,5))
X_test.plot(kind='line', x='ds', y='y', color='blue', ax=ax, label="actual")
X_test_forecast.plot(kind='line', x='ds', y='yhat', color='red', ax=ax, label="forecast")
plt.title("Actual vs prediction (AQI)")
plt.grid(True)
plt.show()<jupyter_output><empty_output><jupyter_text>With this graph we can see that the pred vs actual is some what in the trend, but the range is not being captured ie. the peaks and troughs.<jupyter_code>#collapse_show
#calculating the error

f"The MAPE (Mean absolute Percentage error) for simple Prophet model is: {mean_absolute_per_error(X_test['y'], X_test_forecast['yhat'])}"<jupyter_output><empty_output><jupyter_text>The error % is quite large. We can now optimize the same using hyperparamter tuning.  Different hyperparameters which can be used in the model to optmize the model performance.  - **Holidays** We can add holidays in the model using holiday liberary. Certain custom holidays can also be added.- **n_changepoints** is the number of change happen in the data. Prophet model detects them by its own. By default, its value is 25, which are uniformly placed in the first 80% of the time series. Changing n_changepoints can add value to the model.- **changepoint_prior_scale** to indicate how flexible the changepoints are allowed to be. In other words, how much can the changepoints fit to the data. If you make it high it will be more flexible, but you can end up overfitting. By default, this parameter is set to 0.05- **seasonality_mode** There are 2 types model seasonality mode. Additive & multiplicaticative. By default Prophet fits additive seasonalities, meaning the effect of the seasonality is added to the trend to get the forecast. Prophet can model multiplicative seasonality by setting seasonality_mode='multiplicative' in the model.- **holiday_prior_scale** just like changepoint_prior_scale, holiday_prior_scale is used to smoothning the effect of holidays. By default its value is 10, which provides very little regularization. Reducing this parameter dampens holiday effects- **Seasonalities with fourier_order** Prophet model, by default finds the seasonalities and adds the default parameters of the seasonality. We can modify the seasonalities effect by adding custom seasonalities as add_seasonality in the model with different fourier order.By default Prophet uses a Fourier order of 3 for weekly seasonality and 10 for yearly seasonality.  Hyperparameter tuning using Parameter Grid<jupyter_code>#collapse_show
from sklearn.model_selection import ParameterGrid
params_grid={'seasonality_mode':('multiplicative','additive'),
            'changepoint_prior_scale':[0.1,0.2,0.3,0.4,0.5],
            'n_changepoints':[100,150,200]}
grid=ParameterGrid(params_grid)
cnt=0
for p in grid:
    cnt+=1
f"the tolal number of possible model combinations are: {cnt}"<jupyter_output><empty_output><jupyter_text>Model Tuning<jupyter_code>#collapse_show
start=pd.to_datetime('2019-12-01')
end=pd.to_datetime('2020-07-01')
model_parameters= pd.DataFrame(columns=['MAPE','Parameters'])
for p in grid:
    test= pd.DataFrame()
    print(p)
    train_model=Prophet(changepoint_prior_scale=p['changepoint_prior_scale'],
                       n_changepoints=p['n_changepoints'],
                       seasonality_mode=p['seasonality_mode'],
                       weekly_seasonality=True,
                       daily_seasonality=True,
                       yearly_seasonality= True, interval_width=0.95)
    train_model.fit(X_train)
    train_forecast= train_model.make_future_dataframe(periods=213, freq='D', include_history=False)
    train_forecast=train_model.predict(train_forecast)
    test=train_forecast[['ds','yhat']]
    Actual=dfd[(dfd['ds']>=start) & (dfd['ds']<end)]
    # test=test[(test['ds']>=start) | (test['ds']<end)]
    MAPE=mean_absolute_per_error(Actual['y'],abs(test['yhat']))
    print(' Mean absolute Percentage error-------------------', MAPE)
    model_parameters=model_parameters.append({'MAPE':MAPE, 'Parameters':p}, ignore_index=True)
    
#collapse_show
# Filtering the best parameters based on MAPE

model_parameters.sort_values(by='MAPE').head()
#collapse_show
f"The best parameters are: {model_parameters.loc[0]['Parameters']}"
#collapse_show
#Using the best parameter for model creation 

final_model=Prophet(changepoint_prior_scale=0.1, 
                   n_changepoints=100, 
                   seasonality_mode='multiplicative',
                   weekly_seasonality=True,
                   daily_seasonality= True,
                   yearly_seasonality= True,
                   interval_width=0.95)
#collapse_show
final_model.fit(X_train)
#collapse_show
future=final_model.make_future_dataframe(periods=602, freq='D')
forecast=final_model.predict(future)
forecast[['ds','yhat','yhat_lower','yhat_upper']].tail(10)
#collapse_show
#plotting the model components

fig= final_model.plot_components(forecast)
#collapse_show
# predicting the test values with this model

X_test_final= final_model.predict(X_test)
X_test_final.tail()
#collapse_show
#calculating the error

f"The MAPE (Mean absolute Percentage error) for final Prophet model is: {mean_absolute_per_error(X_test['y'], X_test_final['yhat'])}"<jupyter_output><empty_output><jupyter_text>The error on the test set has gone down, this indicates better data fittment<jupyter_code>#collapse_show
#mapping the test and forcast data

f,ax= plt.subplots(figsize=(12,5))
X_test.plot(kind='line', x='ds', y='y', label="Actual", ax=ax)
X_test_final.plot(kind='line', x='ds', y='yhat', label="forecast", ax=ax)
plt.grid(True)
plt.show()<jupyter_output><empty_output><jupyter_text>Time series decomposition There are three componenets of time series.- Trend: Consistant up/down slope of a time series- Seasonality: Periodic pattern of the time series- Noise: Outliers or missing values<jupyter_code>#collapse_show
# Understanding these factors via plotting graph
X_train.plot(kind='line', x='ds', y='y')
plt.title("Plot of AQI")
plt.grid(True)
#collapse_show
#breaking the plot in different sub components

import statsmodels.api as sm
decomposed_AQI=sm.tsa.seasonal_decompose(X_train['y'], freq=360)
decomposed_AQI.plot()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12268787996378218, 'fraction_numerical': 0.012676238520243176, 'mean_word_length': 6.106158088235294, 'pattern_counts': {'":': 0, '<': 96, '<?xml version=': 0, '>': 97, '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': '23', 'dir': 'jupyter-structured-clean-dedup', 'id': '532333', 'license': 'Apache-2.0', 'n_tokens_mistral': 5368, 'n_tokens_neox': 5055, 'n_words': 1405, 'path': '_notebooks/2020-12-02-Time_Series_analysis_delhi.ipynb', 'repo_name': 'kaushalbundel/blog', 'size': '15462'}