ayashabbar / covid19-ts-model Goto Github PK

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from statsmodels.tsa.vector_ar.var_model import VAR
from statsmodels.tsa.stattools import adfuller

death_data = pd.read_csv('deaths.csv')
confirmed_data = pd.read_csv('confirmed.csv')
death_data = death_data.drop('Province/State', axis = 1)
confirmed_data = confirmed_data.drop('Province/State', axis = 1)
location_data = pd.DataFrame({'Country': confirmed_data['Country/Region'], 'Lat': confirmed_data['Lat'], 'Long' : confirmed_data['Long']}, columns=['Country', 'Lat', 'Long'])
location_data = location_data.groupby('Country').agg('mean')

death_data = death_data.drop('Lat', axis = 1)
death_data = death_data.drop('Long', axis = 1)
confirmed_data = confirmed_data.drop('Lat', axis = 1)
confirmed_data = confirmed_data.drop('Long', axis = 1)

confirmed_data = confirmed_data.groupby('Country/Region').agg('sum')
#confirmed_data = confirmed_data.join(location_data)
death_data = death_data.groupby('Country/Region').agg('sum')
#death_data = death_data.join(location_data)

globalnum1 = confirmed_data.sum(axis=0).rename("Global")
globalnum2 = death_data.sum(axis=0).rename("Global")
countryfilter = confirmed_data.index.isin(["South Korea", "Italy", "China", "US", "Japan", "United Kingdom", "France", "Korea South", "Global"])
confirmed_data = confirmed_data[countryfilter]
death_data = death_data[countryfilter]

# Dividie by population/100000 -> Each point represents 1/1000 of a percent
China = confirmed_data.iloc[0].divide(14377.22468)
France = confirmed_data.iloc[1].divide(652.73511)[2:]
Italy = confirmed_data.iloc[2].divide(604.87167)[9:]
Japan = confirmed_data.iloc[3].divide(1265.87521)
USA = confirmed_data.iloc[4].divide(3304.41184)
Korea = confirmed_data.iloc[5].divide(512.56069)

confirmed_data = confirmed_data.append(globalnum1)
confirmed_data.head()

death_data = death_data.append(globalnum2)
death_data.head()

fig,a =  plt.subplots(2,3)
x = np.arange(1,5)
fig.tight_layout()
a[0][0].plot(range(1,56), confirmed_data.iloc[0])
a[0][0].set_title('China')
a[0][1].plot(range(1,56), confirmed_data.iloc[1])
a[0][1].set_title('France')
a[1][0].plot(range(1,56), confirmed_data.iloc[2])
a[1][0].set_title('Italy')
a[1][1].plot(range(1,56), confirmed_data.iloc[3])
a[1][1].set_title('Japan')
a[0][2].plot(range(1,56), confirmed_data.iloc[4])
a[0][2].set_title('USA')
a[1][2].plot(range(1,56), confirmed_data.iloc[5])
a[1][2].set_title('Korea')
plt.show()

# Dickey Fuller Test for stationarity
def test_stationarity(timeseries):
    
    #Determing rolling statistics
    ema = timeseries.ewm(4).mean()
    rolmean = timeseries.rolling(7).mean()
    stat = rolmean-ema
    rolstd = timeseries.rolling(7).std()
    #Plot rolling statistics:
    plt.plot(timeseries, color='blue',label='Original')
    #plt.plot(rolmean, color='red', label='Rolling Mean')
    plt.plot(stat, color= 'green', label = 'EXP Smoothed')
    plt.plot(rolstd, color='black', label = 'Rolling Std')
    plt.legend(loc='best')
    plt.title('Rolling Mean & Standard Deviation')
    plt.show()
    #Perform Dickey-Fuller test:
    print ('Results of Dickey-Fuller Test:')
    dftest = adfuller(timeseries, 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
    return ema

# Make timeseries stationary
def make_stationary(timeseries):
    new_timeseries = []
    ema = timeseries.ewm(4).mean()
    rolmean = timeseries.rolling(7).mean()
    stat = rolmean-ema
    for i in range(0, 6):
        stat[i] = timeseries[i]
    return stat

def get_ema(timeseries):
    return timeseries.ewm(4).mean()

test_stationarity(Japan)

ewi = pd.read_csv('EWI.csv').transpose()
spy = pd.read_csv('SPY.csv').transpose()
ewq = pd.read_csv('EWQ.csv').transpose()
mchi = pd.read_csv('MCHI.csv').transpose()
ewj = pd.read_csv('EWJ.csv').transpose()

cols = range(0,214)
ewi.drop(ewi.columns[cols],axis=1,inplace=True)
spy.drop(spy.columns[cols],axis=1,inplace=True)
ewq.drop(ewq.columns[cols],axis=1,inplace=True)
mchi.drop(mchi.columns[cols],axis=1,inplace=True)
ewj.drop(ewj.columns[cols],axis=1,inplace=True)
ewj.drop([ewj.index[1], ewj.index[2], ewj.index[3], ewj.index[4]], axis = 0, inplace=True)
spy.drop([spy.index[1], spy.index[2], spy.index[3], spy.index[4]], axis = 0, inplace=True)
ewi.drop([ewi.index[1], ewi.index[2], ewi.index[3], ewi.index[4]], axis = 0, inplace=True)
mchi.drop([mchi.index[1], mchi.index[2], mchi.index[3], mchi.index[4]], axis = 0, inplace=True)
ewq.drop([ewq.index[1], ewq.index[2], ewq.index[3], ewq.index[4]], axis = 0, inplace=True)

ewi = ewi.replace('\d\d\d\d-', '', regex=True).replace('\d1-', '1-', regex=True).replace('\d2-', '2-', regex=True).replace('\d3-', '3-', regex=True).replace('-', '/', regex=True)
ewi.iloc[0] = ewi.iloc[0]+'/20'
mchi = mchi.replace('\d\d\d\d-', '', regex=True).replace('\d1-', '1-', regex=True).replace('\d2-', '2-', regex=True).replace('\d3-', '3-', regex=True).replace('-', '/', regex=True)
mchi.iloc[0] = mchi.iloc[0]+'/20'
ewj = ewj.replace('\d\d\d\d-', '', regex=True).replace('\d1-', '1-', regex=True).replace('\d2-', '2-', regex=True).replace('\d3-', '3-', regex=True).replace('-', '/', regex=True)
ewj.iloc[0] = ewj.iloc[0]+'/20'
spy = spy.replace('\d\d\d\d-', '', regex=True).replace('\d1-', '1-', regex=True).replace('\d2-', '2-', regex=True).replace('\d3-', '3-', regex=True).replace('-', '/', regex=True)
spy.iloc[0] = spy.iloc[0]+'/20'
ewq = ewq.replace('\d\d\d\d-', '', regex=True).replace('\d1-', '1-', regex=True).replace('\d2-', '2-', regex=True).replace('\d3-', '3-', regex=True).replace('-', '/', regex=True)
ewq.iloc[0] = ewq.iloc[0]+'/20'

ewi.columns = ewi.iloc[0]
ewi = ewi.drop('Date',axis = 0)
ewi = ewi.drop(['3/17/20','3/18/20'], axis = 1)

ewq.columns = ewq.iloc[0]
ewq = ewq.drop('Date',axis = 0)
ewq = ewq.drop(['3/17/20','3/18/20'], axis = 1)

ewj.columns = ewj.iloc[0]
ewj = ewj.drop('Date',axis = 0)
ewj = ewj.drop(['3/17/20','3/18/20'], axis = 1)

spy.columns = spy.iloc[0]
spy = spy.drop('Date',axis = 0)
spy = spy.drop(['3/17/20','3/18/20'], axis = 1)

mchi.columns = mchi.iloc[0]
mchi = mchi.drop('Date',axis = 0)
mchi = mchi.drop(['3/17/20','3/18/20'], axis = 1)

fig,a =  plt.subplots(2,3)
x = np.arange(1,5)
fig.tight_layout()
a[0][0].plot(range(1,39), mchi.iloc[0])
a[0][0].set_title('China')
a[0][1].plot(range(1,39), ewq.iloc[0])
a[0][1].set_title('France')
a[1][0].plot(range(1,39), ewi.iloc[0])
a[1][0].set_title('Italy')
a[1][1].plot(range(1,39), ewj.iloc[0])
a[1][1].set_title('Japan')
a[0][2].plot(range(1,39), spy.iloc[0])
a[0][2].set_title('USA')
plt.show()

confirmed_data.transpose().columns.name='Date'

def find_increment(prev, next, count):
    return ((next-prev)/(count+1))

def fill_series(series):
    prev = 0;
    next = series[0];
    count = 0;
    i = 0;
    while i < series.size:
        if series[i] == 0:
            count = count + 1
        elif count != 0:
            prev = next
            next = series[i]
            next = series[i]
            increment = find_increment(prev, next, count)
            add = increment
            while count != 0:
                # set previous values 
                series[i-count] = prev + add
                add += increment
                count = count - 1
            count = 0
            prev = next
            next = series[i]
        else:
            prev = next;
            next = series[i]
        i = i + 1

dates = pd.date_range('01/22/2020', '03/16/2020')
#SPY
spy_price = pd.Series(spy.iloc[0])
spy_price.index = pd.DatetimeIndex(spy_price.index)
spy_price = spy_price.reindex(dates, fill_value=0)
fill_series(spy_price)
spy_vol = pd.Series(spy.iloc[1])
spy_vol.index = pd.DatetimeIndex(spy_vol.index)
spy_vol = spy_vol.reindex(dates, fill_value=0)
fill_series(spy_price)
fill_series(spy_vol)
#EWQ
ewq_price = pd.Series(ewq.iloc[0])
ewq_price.index = pd.DatetimeIndex(ewq_price.index)
ewq_price = ewq_price.reindex(dates, fill_value=0)
fill_series(ewq_price)
ewq_vol = pd.Series(ewq.iloc[1])
ewq_vol.index = pd.DatetimeIndex(ewq_vol.index)
ewq_vol = ewq_vol.reindex(dates, fill_value=0)
fill_series(ewq_vol)
#EWJ
ewj_price = pd.Series(ewj.iloc[0])
ewj_price.index = pd.DatetimeIndex(ewj_price.index)
ewj_price = ewj_price.reindex(dates, fill_value=0)
fill_series(ewj_price)
ewj_vol = pd.Series(ewj.iloc[1])
ewj_vol.index = pd.DatetimeIndex(ewj_vol.index)
ewj_vol = ewj_vol.reindex(dates, fill_value=0)
fill_series(ewj_vol)
#EWI
ewi_price = pd.Series(ewi.iloc[0])
ewi_price.index = pd.DatetimeIndex(ewi_price.index)
ewi_price = ewi_price.reindex(dates, fill_value=0)
fill_series(ewi_price)
ewi_vol = pd.Series(ewi.iloc[1])
ewi_vol.index = pd.DatetimeIndex(ewi_vol.index)
ewi_vol = ewi_vol.reindex(dates, fill_value=0)
fill_series(ewi_vol)
#MCHI
mchi_price = pd.Series(mchi.iloc[0])
mchi_price.index = pd.DatetimeIndex(mchi_price.index)
mchi_price = mchi_price.reindex(dates, fill_value=0)
fill_series(mchi_price)
mchi_vol = pd.Series(mchi.iloc[1])
mchi_vol.index = pd.DatetimeIndex(mchi_vol.index)
mchi_vol = mchi_vol.reindex(dates, fill_value=0)
fill_series(mchi_vol)

spy = pd.concat([spy_price, spy_vol], axis=1)
ewq = pd.concat([ewq_price, ewq_vol], axis=1)
ewi = pd.concat([ewi_price, ewi_vol], axis=1)
ewj = pd.concat([ewj_price, ewj_vol], axis=1)
mchi = pd.concat([mchi_price, mchi_vol], axis=1)

# Reasonably stationary for EWJ, but not too stationary for SPY and MCHI
# Stretch goal: Try other forms of smoothing; exponential smoothing generates too many sharp changes in SPY/MCHI
test_stationarity(spy.iloc[:,1])
test_stationarity(spy.iloc[:,0])
test_stationarity(ewj.iloc[:,1])
test_stationarity(ewj.iloc[:,0])
test_stationarity(mchi.iloc[:,1])
test_stationarity(mchi.iloc[:,0])

# Modeling Japan
data = pd.concat([make_stationary(ewj.iloc[:,0]), make_stationary(ewj.iloc[:,1]), make_stationary(death_data.iloc[3]),
                  make_stationary(Japan).astype(float)], axis=1)

# Fitting the VAR Model
model = VAR(data)
model_fit = model.fit()
ewj_pred = model_fit.forecast(model_fit.y, steps=1)

ewj_pred[0][0] + get_ema(ewj.iloc[:,0])[-1] # Predicted 51.1304 EWJ

corona_data = confirmed_data.transpose()
df = data = pd.concat([ewj.iloc[:,0].astype(float), ewi.iloc[:,0].astype(float), ewq.iloc[:,0].astype(float),
                       mchi.iloc[:,0].astype(float), spy.iloc[:,0].astype(float), 
                      corona_data.iloc[:,0], corona_data.iloc[:,1], corona_data.iloc[:,2], corona_data.iloc[:,3], 
                      corona_data.iloc[:,4], death_data.transpose().iloc[:,6]], axis=1)
#df = df.rename(columns ={"Global" : "Deaths"})
#df.rename(columns={df.columns[0]: "EWJ" }, inplace = True)
df.columns = ['EWJ', 'EWI','EWQ','MCHI', 'SPY', 'China', 'France', 'Italy', 'Japan', 'USA', 'Global Deaths']

df

corr = df.corr()
fig, ax = plt.subplots(figsize=(10,10))   
sns.heatmap(corr, annot = True, linewidths = 0.1, ax = ax)