dsc_project.py

# -*- coding: utf-8 -*-
"""DSc_Project.ipynb

Automatically generated by Colaboratory.

Original file is located at
    https://colab.research.google.com/drive/1f4EFpo3-vdkiqX56rapp3if5TkcSEqse
"""

from google.colab import drive
drive.mount('/content/drive')

# Commented out IPython magic to ensure Python compatibility.
import numpy as np 
import pandas as pd 
import seaborn as sns 
from fbprophet import *
import itertools
from matplotlib import pyplot as plt 
# %matplotlib inline
sns.set_style("whitegrid")
import statsmodels.api as sm
from sklearn.linear_model import LogisticRegression,LinearRegression
from sklearn.ensemble import RandomForestClassifier,RandomForestRegressor
from sklearn.neighbors import KNeighborsClassifier,KNeighborsRegressor
from sklearn.tree import DecisionTreeClassifier
from sklearn.svm import SVC, LinearSVC
import string

from sklearn.metrics import *
from sklearn.model_selection import GridSearchCV
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import LabelEncoder

import numpy as np # linear algebra
import pandas as pd # pandas for dataframe based data processing and CSV file I/O


import statistics as stats
from pandas.plotting import lag_plot
from pandas.plotting import autocorrelation_plot

from statsmodels.tsa.arima_model import ARIMA
import statsmodels.api as sm
import warnings
warnings.filterwarnings('ignore')

# df=pd.read_csv('data.csv',encoding = "ISO-8859-1")
df=pd.read_csv('drive/My Drive/DSc/data.csv',encoding = "ISO-8859-1")

df.drop(['stn_code','agency','sampling_date','location_monitoring_station'],axis=1,inplace=True)
data=df
data['date'] = pd.to_datetime(data['date'],format='%Y-%m-%d') # date parse
data['year'] = data['date'].dt.year # year
# data['year'] = data['year'].fillna(0.0).astype(int)
data = data[(data['year']>0)]
data.info()

# from sklearn.preprocessing import Imputer
# imputer = Imputer(missing_values = 'NaN', strategy = 'mean', axis = 0)
# imputer = imputer.fit(data.iloc[:, 3:8].values)
# data.iloc[:,3:8] = imputer.transform(data.iloc[:, 3:8].values)
# data = data.apply (pd.to_numeric, errors='coerce')
data=data.drop(columns=['pm2_5','rspm','spm'])
data = data.dropna()
data.info()
data.head()

data['type'].describe()
#With 10 Unique labels, we will fill the null values by the most common type, which is 'Residential, Rural and Other Areas'.
common_value='Residential,Rural and other Areas'
data['type']=data['type'].fillna(common_value)
data.info()

pd.isnull(data).sum()

f, ax = plt.subplots(figsize=(15,15))
ax.set_title('{} by state and year'.format('so2'))
sns.heatmap(data.pivot_table('so2', index='state',
                columns=['year'],aggfunc='median',margins=True),
                annot=True,cmap="BuPu", linewidths=.5, ax=ax,cbar_kws={'label': 'Annual Average'})

def high_lvl(indicator="so2"):
    plt.figure(figsize=(10,5))
    ind = data[[indicator, 'location', 'state', 'date']].groupby('state', as_index=False).max()
    highest = sns.barplot(x='state', y=indicator, data=ind)
    highest.set_title("Highest ever {} levels recorded by state".format(indicator))
    plt.xticks(rotation=90)
high_lvl()

def mean_lvl(indicator="so2"):
    plt.figure(figsize=(10,5))
    ind = data[[indicator, 'location', 'state', 'date']].groupby('state', as_index=False).mean()
    highest = sns.barplot(x='state', y=indicator, data=ind)
    highest.set_title("Mean ever {} levels recorded by state".format(indicator))
    plt.xticks(rotation=90)
mean_lvl()

data.keys()

def pro(d1,x,m):
  d1['date']=pd.to_datetime(d1['date'])
  d1=d1.rename(columns={"date": "ds", x: "y"})
  m.fit(d1)
  future = m.make_future_dataframe(periods=365*4)
  future.tail()
  forecast = m.predict(future)
  forecast[['ds', 'yhat', 'yhat_lower', 'yhat_upper']].tail()
  fig1 = m.plot(forecast)
  fig2 = m.plot_components(forecast)

m=Prophet()
dso2=data.drop(['state', 'location', 'type', 'no2', 'rspm', 'spm', 'pm2_5',
        'year'],axis=1)
pro(dso2,"so2",m)

m=Prophet()
dno2=data.drop(['state', 'location', 'type', 'so2', 'rspm', 'spm', 'pm2_5',
        'year'],axis=1)
pro(dno2,"no2",m)

m=Prophet()

dno2delhi=data[data['location']=='Delhi']
# print(dno2delhi)
# dno2delhi=dno2delhi.drop(['state', 'location','type', 'so2', 'rspm', 'spm', 'pm2_5','year'],axis=1)
dno2delhi=dno2delhi.drop(['state', 'location','type', 'so2', 'year'],axis=1)

pro(dno2delhi,"no2",m)

# future = m.make_future_dataframe(periods=365)
# forecast = m.predict(future)
# m.plot_components(forecast)

m=Prophet()

dno2mum=data[data['location']=='Mumbai']
# print(dno2delhi)
# dno2mum=dno2mum.drop(['state', 'location','type', 'so2', 'rspm', 'spm', 'pm2_5','year'],axis=1)
dno2mum=dno2mum.drop(['state', 'location','type', 'so2', 'year'],axis=1)
pro(dno2mum,"no2",m)

pip install pyflux

df = data[['no2','year','state']].groupby(["year"]).median().reset_index().sort_values(by='year',ascending=False)
f,ax=plt.subplots(figsize=(15,5))
sns.pointplot(x='year', y='no2', data=df).set_title("NO2 Analysis")

dk = data[['so2','year','state']].groupby(["year"]).median().reset_index().sort_values(by='year',ascending=False)
f,ax=plt.subplots(figsize=(15,5))
sns.pointplot(x='year', y='so2', data=dk).set_title("SO2 Analysis")

print(dno2delhi.info() )
autocorrelation_plot(dno2delhi['no2'])

m=Prophet()

dso2delhi=data[data['location']=='Delhi']
# print(dno2delhi)
# dso2delhi=dso2delhi.drop(['state', 'location','type', 'no2', 'rspm', 'spm', 'pm2_5','year'],axis=1)
dso2delhi=dso2delhi.drop(['state', 'location','type', 'no2','year'],axis=1)
pro(dso2delhi,"so2",m)

X = dk['so2'].values
size = int(len(X) * 0.60)
train, test = X[0:size], X[size:len(X)]
history = [x for x in train]
predictions = list()
for t in range(len(test)):
    model = ARIMA(history, order=(5,1,0))
    model_fit = model.fit(disp=0)
    obs = test[t]
    output = model_fit.forecast()
    yhat = output[0]
    predictions.append(yhat)
    history.append(obs)
    print('predicted=%f, expected=%f' % (yhat, obs))
error = np.sqrt(mean_squared_error(test, predictions))
print('Test RMSE: %.5f' % error)
# plot
# plt.plot(train)
plt.plot(test)
plt.plot(predictions, color='red')
plt.show()

# print(df.info())
X = df['no2'].values
dat=df['year'].values
size = int(len(X) * 0.60)
train, test = X[0:size], X[size:len(X)]
history = [x for x in train]
predictions = list()
for t in range(len(test)):
    model = ARIMA(history, order=(1,1,0))
    model_fit = model.fit(disp=0,maxiter=1000)
    obs = test[t]
    output = model_fit.forecast()
    yhat = output[0]
    predictions.append(yhat)
    history.append(obs)
    print('predicted=%f, expected=%f' % (yhat, obs))
error = np.sqrt(mean_squared_error(test, predictions))
print('Test RMSE: %.5f' % error)
# plot
# plt.plot(dat[0:size],train)
plt.plot(dat[size:len(X)],test,color='green')
plt.plot(dat[size:len(X)],predictions, color='red')
plt.show()
print(model_fit.summary2())

X_diff = X - df['no2'].shift()
model = ARIMA(X, order=(3, 1, 1))
results_MA = model.fit(disp=-1)  
plt.plot(X_diff)
plt.plot(results_MA.fittedvalues, color='red')
print(results_MA.summary2())

p = d = q = range(0, 5)
pdq = list(itertools.product(p, d, q))
seasonal_pdq = [(x[0], x[1], x[2], 12) for x in list(itertools.product(p, d, q))]
print('Examples of parameter combinations for Seasonal ARIMA...')
print('SARIMAX: {} x {}'.format(pdq[1], seasonal_pdq[1]))
print('SARIMAX: {} x {}'.format(pdq[1], seasonal_pdq[2]))
print('SARIMAX: {} x {}'.format(pdq[2], seasonal_pdq[3]))

print('SARIMAX: {} x {}'.format(pdq[2], seasonal_pdq[4]))

lis=[]
for param in pdq:
    for param_seasonal in seasonal_pdq:
        try:
            mod = sm.tsa.statespace.SARIMAX(X_diff,order=param,seasonal_order=param_seasonal,enforce_stationarity=False)
            results = mod.fit()
            print('SARIMAX{}x{}12 - AIC:{}'.format(param, param_seasonal, results.aic))
            lis.append( results.aic)
        except:
            continue
# lis.sort()
# print(lis)

np.sort(lis)[0:40]

# ARIMA(2, 0, 1)x(2, 0, 0, 12)12
mod = sm.tsa.statespace.SARIMAX(X_diff,
                                order=(2, 0, 1),
                                seasonal_order=(2, 0, 0, 12),
                                enforce_stationarity=False,
                                enforce_invertibility=False)
results = mod.fit()
print(results.summary().tables[1])

results.plot_diagnostics()
plt.show()

def top_and_bottom_10_states(indicator="so2"):
    fig, ax = plt.subplots(2,1, figsize=(20, 12))
    
    ind = data[[indicator, 'state']].groupby('state', as_index=False).mean().sort_values(by=indicator,ascending=False)
    
    top10 = sns.barplot(x='state', y=indicator, data=ind[:10], ax=ax[0], color='R')
    top10.set_title("Top 10 states by {} (1991-2016)".format(indicator))
    top10.set_ylabel("so2 (µg/m3)")
    top10.set_xlabel("State")
    
    bottom10 = sns.barplot(x='state', y=indicator, data=ind[-10:], ax=ax[1], color='g')
    bottom10.set_title("Bottom 10 states by {} (1991-2016)".format(indicator))
    bottom10.set_ylabel("so2 (µg/m3)")
    bottom10.set_xlabel("State")

top_and_bottom_10_states("no2")

top_and_bottom_10_states("so2")

itop_and_bottom_10_states("rspm")

"""End

---



"""