# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# + id="uCcCj6_oMOH2" colab_type="code" outputId="5fcdee34-e4e8-4439-b59d-f14c2a0eb649" executionInfo={"status": "ok", "timestamp": 1589311073939, "user_tz": 240, "elapsed": 18561, "user": {"displayName": "<NAME>", "photoUrl": "", "userId": "16347493411284187623"}} colab={"base_uri": "https://localhost:8080/", "height": 122}
from google.colab import drive
drive.mount('/content/gdrive')

# + id="oppbm5AkMak0" colab_type="code" outputId="dd3ccd16-c71f-47be-de90-788694b2501b" executionInfo={"status": "ok", "timestamp": 1589311212948, "user_tz": 240, "elapsed": 893, "user": {"displayName": "<NAME>", "photoUrl": "", "userId": "16347493411284187623"}} colab={"base_uri": "https://localhost:8080/", "height": 34}
# %cd '/content/gdrive/Shared drives/Cse 544'

# + id="WCSxptVMMeLT" colab_type="code" outputId="acc1f2d2-c4a0-410e-ac75-930af22d3141" executionInfo={"status": "ok", "timestamp": 1589311215167, "user_tz": 240, "elapsed": 1674, "user": {"displayName": "<NAME>", "photoUrl": "", "userId": "16347493411284187623"}} colab={"base_uri": "https://localhost:8080/", "height": 34}
# !ls

# + id="wLQNVPyHL8gn" colab_type="code" colab={}
import pandas as pd
import matplotlib.pyplot as plt
from datetime import datetime, timedelta
import numpy as np
import math
from pprint import pprint
from scipy.stats import geom, poisson, binom
from scipy.stats import chi2_contingency

states = ['California', 'New York', 'Washington']


# + id="B5n5XBiSL8gr" colab_type="code" colab={}
def outlier(data, length, alpha):
    quartile1 = data[length]
    quartile2 = data[length*2]
    quartile3 = data[length*3]
    interquartile = quartile3 - quartile1
    lower_limit = quartile1 - (alpha * interquartile)
    upper_limit = quartile3 + (alpha * interquartile)
    print("Quartile1: " + str(quartile1))
    print("Quartile2: " + str(quartile2))
    print("Quartile3: " + str(quartile3))
    print("Interquartile: " + str(interquartile))
    print("Lower limit: " + str(lower_limit))
    print("Upper limit: " + str(upper_limit))
    flag = 0
    for i,num in enumerate(data):
        if num < lower_limit or num > upper_limit: 
            flag = +1
            print("Excluded value: " + str(num))
    return flag, lower_limit, upper_limit


# + id="F_rlSffSL8gt" colab_type="code" colab={}
def retrieve_data(covid_fpath, stock_fpath):
    # Data cleaning and outlier for Covid dataset
    data_states = pd.read_csv(covid_fpath)
    states = ['California', 'New York', 'Washington']
    data_states['date'] = pd.to_datetime(data_states["date"])

    # Generating all dates to find missing values
    all_dates = set()
    start_date = datetime.strptime('2020-01-01', '%Y-%m-%d')
    end_date = data_states.date.max()
    while start_date<end_date:
        all_dates.add(start_date)
        start_date = start_date + timedelta(1)

    state_df = pd.DataFrame(columns=['date', 'state', 'cases', 'deaths'])

    for st in states:
        st_df = data_states[data_states.state == st][['date', 'state', 'cases', 'deaths']].sort_values(by='date', axis=0, ascending=True)
        st_missing_dates = all_dates - set(st_df.date.tolist())
        new_df = []
        # Replace missing values with 0 (cases, deaths)
        for dt in st_missing_dates:
            new_df.append([dt, st, 0, 0])
        new_df = pd.DataFrame(new_df, columns=['date', 'state', 'cases', 'deaths'])

        # Converting cumulative case counts to new cases
        st_new_cases = np.array([0]*st_df.shape[0])
        old_case_count = 0 
        for i in range(st_df.shape[0]):
            st_new_cases[i] = st_df.iloc[i].cases - old_case_count
            old_case_count = st_df.iloc[i].cases
        st_df['cases'] = st_new_cases

        # Converting cumulative death counts to new deaths
        st_new_deaths = np.array([0]*st_df.shape[0])
        old_death_count = 0
        for i in range(st_df.shape[0]):
            st_new_deaths[i] = st_df.iloc[i].deaths - old_death_count
            old_death_count = st_df.iloc[i].deaths
        st_df['deaths'] = st_new_deaths
        state_df = state_df.append(new_df)
        state_df = state_df.append(st_df)

    # contains data of 3 states (at date, state level)
    state_df.sort_values(by=['date', 'state'], inplace=True)
    state_df.reset_index(drop=True, inplace=True)

    # contains daily data of 3 states (at date level)
    current_df = state_df[['date', 'cases', 'deaths']]
    current_df = current_df.groupby("date").sum()
    current_df['date'] = current_df.index
    current_df.reset_index(drop=True, inplace=True)
    current_df = current_df[['date', 'cases', 'deaths']]
    
    # alpha for outlier detection
    outlier_alpha = 1.4
    print("==========Cases Outlier===========")
    domain = current_df['cases'].sort_values().values.tolist()
    case_flag, case_lower, case_upper = outlier(domain, len(domain)//4, outlier_alpha)
    print("Outliers in cases: ", case_flag)
    
    print("==========Deaths Outlier===========")
    domain = current_df['deaths'].sort_values().values.tolist()
    print(current_df[current_df['deaths'] == 892])
    death_flag, death_lower, death_upper = outlier(domain, len(domain)//4, outlier_alpha)
    if death_flag == 0: print("No Outlier detected")
    print("Outliers in deaths: ", death_flag)
    
    # Removing outliers from covid dataset
    current_df = current_df.loc[(current_df['cases']>=case_lower) & (current_df['cases']<=case_upper) \
                  & (current_df['deaths']>=death_lower) & (current_df['deaths']<=death_upper)].reset_index(drop=True)
    
    # Stock dataset cleaning containing stock prices after 2020-01-01
    stock_df = pd.read_csv(stock_fpath, names=['date', 'Open', 'High', 'Low', 'Close', 'Change'], header=0)
    stock_df = stock_df[stock_df.Open.isna()==False][['date', 'Close', 'Change']]
    stock_df['date'] = pd.to_datetime(stock_df["date"])
    end_date = current_df.date.max()
    stock_df = stock_df.loc[(stock_df.date>='2020-01-01') & (stock_df.date<=end_date)].sort_values(by='date').reset_index(drop=True)
    
    print("==========Stock Outlier===========")
    domain = stock_df['Close'].sort_values().values.tolist()
    stock_flag, stock_lower, stock_upper = outlier(domain, len(domain)//4, outlier_alpha)
    print("Outliers in stock: ", stock_flag)
    
    # Removing outliers from stock dataset
    stock_df = stock_df.loc[(stock_df['Close']>=stock_lower) & (stock_df['Close']<=stock_upper)].reset_index(drop=True)

    # Combined dataset containing both stock and covid in single dataframe (inner join)
    covid_stock_df = current_df.join(stock_df.set_index('date'), on='date', how='inner', lsuffix='covid_', rsuffix='stock_').reset_index(drop=True)
    
    # complete stock dataset
    all_stock_df = pd.read_csv(stock_fpath, names=['date', 'Open', 'High', 'Low', 'Close', 'Change'], header=0)
    all_stock_df = all_stock_df[all_stock_df.Open.isna()==False][['date', 'Close', 'Change']]
    all_stock_df['date'] = pd.to_datetime(all_stock_df["date"])
    all_stock_df = all_stock_df.reset_index(drop=True)
    current_df = current_df[current_df['date'] <= "2020-04-30"]
    stock_df = stock_df[stock_df['date'] <= "2020-04-30"]
    covid_stock_df = covid_stock_df[covid_stock_df['date'] <= "2020-04-30"]
    state_df = state_df[state_df['date'] <= "2020-04-30"]
    all_stock_df = all_stock_df[all_stock_df['date'] <= "2020-04-30"]
    return current_df, stock_df, covid_stock_df, state_df, all_stock_df



# + [markdown] id="EOj2ZrrbSYB8" colab_type="text"
# # Data retrieval and cleaning
#
#

# + id="j1OkTF2BL8gw" colab_type="code" outputId="926c7700-9852-49bf-f453-f41cbd242366" executionInfo={"status": "ok", "timestamp": 1589311801266, "user_tz": 240, "elapsed": 1607, "user": {"displayName": "<NAME>", "photoUrl": "", "userId": "16347493411284187623"}} colab={"base_uri": "https://localhost:8080/", "height": 493}
covid_fpath = 'us-states.csv'
stock_fpath = 'stocks.csv'
current_df, stock_df, covid_stock_df, state_df, all_stock_df = retrieve_data(covid_fpath, stock_fpath)
stock_df.shape, covid_stock_df.shape, all_stock_df.shape

# + [markdown] id="ao0tHLJgqjlx" colab_type="text"
# # Data Visualization

# + id="72YgM5Vyqlf3" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 258} outputId="371ff281-6c2e-4e73-8f84-225fcdb4ccbd" executionInfo={"status": "ok", "timestamp": 1589316844471, "user_tz": 240, "elapsed": 495, "user": {"displayName": "<NAME>", "photoUrl": "", "userId": "16347493411284187623"}}
import matplotlib.pyplot as plt
import seaborn as sns
curr_filter = (covid_stock_df['date']>="2020-03-01") & (covid_stock_df['date']<"2020-05-01")
df = covid_stock_df.loc[curr_filter]
curr_filter = (state_df['date']>="2020-03-01") & (state_df['date']<"2020-05-01")
states = state_df.loc[curr_filter]
df["Market Volatility"] = (round(df['Change']/df['Close'] * 100,2))
df.drop(['Change'], axis=1, inplace=True)

# + [markdown] id="MCEcl-cpqr9v" colab_type="text"
# 1. We filter out the data after March till April End as there were hardly any cases in February for Washnigton and California and there were no cases in New york.

# + id="Ht10xy6tqnIY" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 244} outputId="33f75594-b8d6-4cea-b512-ff915b09c731" executionInfo={"status": "ok", "timestamp": 1589316876281, "user_tz": 240, "elapsed": 1322, "user": {"displayName": "<NAME>", "photoUrl": "", "userId": "16347493411284187623"}}
plt.figure(figsize=(15,3))
for state,color in [('Washington', 'Green'), ('California', 'Red'), ('New York', 'Blue')]:
    x = states.loc[states['state']==state]['date']
    y = states.loc[states['state']==state]['cases']
    plt.plot(x, y, c=color, label=state)
plt.legend()
plt.title("No. of cases")

# + [markdown] id="TIyy5UY3qx_n" colab_type="text"
# The plot above shows the timeline for the number of cases for each state. We can observe the high number of cases for New York as compared to California and Washington.

# + id="4CHplbC1qyxF" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 244} outputId="1d642537-6499-4122-c283-4e942aa2f1b8" executionInfo={"status": "ok", "timestamp": 1589316897603, "user_tz": 240, "elapsed": 1055, "user": {"displayName": "<NAME>", "photoUrl": "", "userId": "16347493411284187623"}}
plt.figure(figsize=(15,3))
for state,color in [('Washington', 'Green'), ('California', 'Red'), ('New York', 'Blue')]:
    x = states.loc[states['state']==state]['date']
    y = states.loc[states['state']==state]['deaths']
    plt.plot(x, y, c=color, label=state)
plt.legend()
# plt.yscale('log')
plt.title("No. of deaths")

# + [markdown] id="sCeosalRq0z_" colab_type="text"
# The plot above shows the timeline for the number of deaths for each state. We can observe the high number of deaths for New York as compared to California and Washington.

# + id="Jt4b9Oznq3IX" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 764} outputId="280e0950-7945-41c3-dcf3-8cb5d2516c1e" executionInfo={"status": "ok", "timestamp": 1589316926728, "user_tz": 240, "elapsed": 4991, "user": {"displayName": "<NAME>", "photoUrl": "", "userId": "16347493411284187623"}}
sns.pairplot(df)


# + [markdown] id="yVqpgAPzq7Bm" colab_type="text"
# **Scatter Plot Matrix**
#
# This kind of matrix is useful to bivariate distributions and univariate(histogram) distribution of a variable easily.
# 1.   From the plot between #cases and #deaths, we see that as number of cases increases, the death count is high as well.
# 1.   From the plot between #cases vs Close Price and #deaths vs Close Price, we can observe that during the initial phase, the stock prices were low, but there seems to be no correlation between #cases and #deaths with the Closing Price.
# 2.   From the above, market volatility seems to be more during Covid.
#
#

# + [markdown] id="VDo13GomjL4J" colab_type="text"
# # Required Inferences 1 (AR and EWMA)

# + [markdown] id="pXt_nY-Xjdzw" colab_type="text"
# **Calculation of $\hat{\beta}$**
#
# This is the code segment for calculating the estimates of the coefficients.i.e. Beta OLS matrix. For AR(3), the X dataset is constructed by taking into account the previous 3 days of observed cases and deaths. Similarly for AR(5), X dataset consiers the previous 5 days of information. The start date of   consideration is April 1, 2020. Thus the first new cases and deaths number observed placed in the Y_obs_cases and Y_obs_deaths is taken for 3 days later i.e. April 4, 2020. After the prearation of the X and Y matrices, the Beta OLS is calculated with the formula: 
# $\hat{\beta} = (X^{T}X)^{-1} X^{T} Y$
#
#
#
#

# + id="tzXF44Q6jPCu" colab_type="code" colab={}
def beta_hat_calc(current_df, iter,n):
  start_date = datetime.strptime('2020-04-01', '%Y-%m-%d')
  start_index = current_df.index[current_df['date'] == start_date][0]
  X_cases = []
  X_deaths = []
  Y_obs_cases=[]
  Y_obs_deaths = []
  col_cases = current_df['cases']
  col_deaths = current_df['deaths']
  if n == 3:
    for i in range(start_index , start_index + iter):
      Y_obs_cases.append(col_cases[i+3])
      Y_obs_deaths.append(col_deaths[i+3])
      X_cases.append([1,col_cases[i+2],col_cases[i+1], col_cases[i]])
      X_deaths.append([1,col_deaths[i+2],col_deaths[i+1], col_deaths[i]])
  if n == 5:
    for i in range(start_index , start_index + iter):
      Y_obs_cases.append(col_cases[i+5])
      Y_obs_deaths.append(col_deaths[i+5])
      X_cases.append([1,col_cases[i+4],col_cases[i+3], col_cases[i+2], col_cases[i+1], col_cases[i]])
      X_deaths.append([1,col_deaths[i+4],col_deaths[i+3], col_deaths[i+2], col_deaths[i+1], col_deaths[i]])
  Y_obs_cases = np.array(Y_obs_cases)
  Y_obs_deaths = np.array(Y_obs_deaths)
  X_cases = np.array(X_cases)
  X_deaths = np.array(X_deaths)
  beta_hat_cases = np.dot(np.dot(np.linalg.pinv(np.dot(X_cases.transpose(),X_cases)),X_cases.transpose()),Y_obs_cases)
  beta_hat_deaths = np.dot(np.dot(np.linalg.pinv(np.dot(X_deaths.transpose(),X_deaths)),X_deaths.transpose()),Y_obs_deaths)
  return (beta_hat_cases, beta_hat_deaths, Y_obs_cases, Y_obs_deaths)


# + [markdown] id="_AaB6NnqjlhW" colab_type="text"
# **MAPE**
#
# The function takes in in the predicted values and their respective error predictions and calculates the mean average percentage error using:
# MAPE% =  $(\sum_{i=0}^{n} |E_i|/ Y_i)/ n $.This MAPE% is then returned.

# + id="YDdR02WXjiZk" colab_type="code" colab={}
def mape(Y_obs,error_hat):
  n=len(error_hat)
  sum=0
  for i in range(n):
    sum+= abs(error_hat[i])/Y_obs[i]
  return (sum*100/n)


# + [markdown] id="vrSz_g2Djt8y" colab_type="text"
# **MSE**
#
# This function takes in the predicted errror values for the 7 days and calculates the mean squared error by applying the formula: 
# $MSE = \sum_{i=1}^{n} (E_i)^2$. This error is then returned for further interpretation.

# + id="GqBFpZXUjwFZ" colab_type="code" colab={}
def mse(error_hat):
  n = len(error_hat)
  sum=0
  for i in range(n):
    sum+= error_hat[i] ** 2
  return (sum/n)


# + [markdown] id="b-jMybV7jzn2" colab_type="text"
# **AR(3)**
#
# The function will make predictions for the i*th* day of the week by taking into consideration all the data from the previous days uptil April 1, 2020. Thus, we will have 17 data points as one is outlier for training for the first day of the fourth week. It calls the `beta_hat_calc()` and uses these values as the coefficients of the previously observed cases and deaths. Thus, $\hat{Y} = \hat{\beta}Y_{obs}$, where the $\hat{\beta}$ has 4 values and $Y_{obs}$ has the values of the three previous days.

# + id="8Q-Tqzyfj2T5" colab_type="code" colab={}
def ar_3(current_df):
  error_hat_cases =[]
  error_hat_deaths = []
  for day in range(7):
    beta_hat_cases , beta_hat_deaths, Y_obs_cases, Y_obs_deaths = beta_hat_calc(current_df,17+day,3)
    Y_pred_cases = beta_hat_cases[0] + (beta_hat_cases[1] * Y_obs_cases[-1]) + (beta_hat_cases[2] * Y_obs_cases[-2]) + (beta_hat_cases[3] * Y_obs_cases[-3])
    Y_pred_deaths = beta_hat_deaths[0] + (beta_hat_deaths[1] * Y_obs_deaths[-1]) + (beta_hat_deaths[2] * Y_obs_deaths[-2]) + (beta_hat_deaths[3] * Y_obs_deaths[-3])
    error_hat_cases.append(Y_obs_cases[-day] - Y_pred_cases)
    error_hat_deaths.append(Y_obs_deaths[-day] - Y_pred_deaths)
  print("MAPE for no. of new cases using AR(3) =" + str(mape(Y_obs_cases[-7:] , error_hat_cases)) + "%")
  print("MSE for no. of new cases using AR(3) = " + str(mse(error_hat_cases)))
  print("MAPE for no. of new deaths using AR(3) =" + str(mape(Y_obs_deaths[-7:] , error_hat_deaths)) + "%")
  print("MSE for no. of new deaths using AR(3) = " + str(mse(error_hat_deaths)))


# + [markdown] id="a5Q24YpLj4aH" colab_type="text"
# **AR(5)**
#
# The function will make predictions for the i*th* day of the week by taking into consideration all the data from the previous days uptil April 1, 2020. Thus, we will have 21 days out of which one is an outlier and thus, 15 data points for training for the first day of the fourth week. It calls the `beta_hat_calc()` and uses these values as the coefficients of the previously observed cases and deaths. Thus, $\hat{Y} = \hat{\beta}Y_{obs}$, where the $\hat{\beta}$ has 6 values and $Y_{obs}$ has the values of the five previous days.

# + id="NVO0uvKdj6Cm" colab_type="code" colab={}
def ar_5(current_df):
  error_hat_cases =[]
  error_hat_deaths = []
  for day in range(7):
    beta_hat_cases , beta_hat_deaths, Y_obs_cases, Y_obs_deaths = beta_hat_calc(current_df,15+day,5)
    Y_pred_cases = beta_hat_cases[0] + (beta_hat_cases[1] * Y_obs_cases[-1]) + (beta_hat_cases[2] * Y_obs_cases[-2]) + (beta_hat_cases[3] * Y_obs_cases[-3]) + (beta_hat_cases[4] * Y_obs_cases[-4]) + (beta_hat_cases[5] * Y_obs_cases[-5])
    Y_pred_deaths = beta_hat_deaths[0] + (beta_hat_deaths[1] * Y_obs_deaths[-1]) + (beta_hat_deaths[2] * Y_obs_deaths[-2]) + (beta_hat_deaths[3] * Y_obs_deaths[-3]) + (beta_hat_deaths[4] * Y_obs_deaths[-4]) + (beta_hat_deaths[5] * Y_obs_deaths[-5])
    error_hat_cases.append(Y_obs_cases[-day] - Y_pred_cases)
    error_hat_deaths.append(Y_obs_deaths[-day] - Y_pred_deaths)
  print("MAPE for no. of new cases using AR(5) =" + str(mape(Y_obs_cases[-7:] , error_hat_cases)) + "%")
  print("MSE for no. of new cases using AR(5) = " + str(mse(error_hat_cases)))
  print("MAPE for no. of new deaths using AR(5) =" + str(mape(Y_obs_deaths[-7:] , error_hat_deaths)) + "%")
  print("MSE for no. of new deaths using AR(5) = " + str(mse(error_hat_deaths)))


# + [markdown] id="cxvfoKIhj8gS" colab_type="text"
# **EWMA**
#
# The Exponentially Weighted Moving Average is a statistic for monitoring the process that averages the data in a way that gives certain weights to data as they are further removed in time. In this function, parameters accepted is the complete data frame and $\alpha$. The calculation for this is done with the formula: $\hat{y}_{t+1|t} = \alpha y_{t} + (1-\alpha)\hat{y}_{t|t-1}$
# When the $\hat{Y}_{4th\_week\_1st\_day}$ data is calculated, the estimated value of $3^{rd}$ week last day is taken aqs the observed value for this day. The MAPE and MSE is then calculated for the values of all 7 days each for no. of cases and deaths. 

# + id="EzEKeGmbj-rE" colab_type="code" colab={}
def ewma(current_df,alpha):
  #declaring the start date which will be the seed value for the ewma 7 day estimation
  start_date = datetime.strptime('2020-04-22', '%Y-%m-%d')
  start_i = current_df.index[current_df['date'] == start_date][0]
  col_cases = current_df['cases'] #the coloumn of new cases in the current_df
  col_deaths = current_df['deaths'] #the coloumn of new deaths in the current_df
  #Initializing the Y lists with the seed value
  Y_obs_cases = [col_cases[start_i]]
  Y_obs_deaths = [col_deaths[start_i]]
  Y_pred_cases = [col_cases[start_i]]
  Y_pred_deaths = [col_deaths[start_i]]
  error_hat_cases = []
  error_hat_deaths = []
  # p_alpha = 1- alpha
  for i in range(7):
    Y_pred_cases.append( (alpha*Y_obs_cases[-1]) + ((1-alpha)*Y_pred_cases[-1]))
    Y_pred_deaths.append( (alpha*Y_obs_deaths[-1]) + ((1-alpha)*Y_pred_deaths[-1]))
    Y_obs_cases.append(col_cases[start_i +i+1])
    Y_obs_deaths.append(col_deaths[start_i +i+1])
    error_hat_cases.append(Y_obs_cases[i+1] - Y_pred_cases[i+1])
    error_hat_deaths.append(Y_obs_deaths[i+1] - Y_pred_deaths[i+1])
  # print(alpha)
  print("MAPE for EWMA for new cases with alpha=" + str(alpha)+ " is " + str(mape(Y_obs_cases[1:] , error_hat_cases)) + "%")
  print("MSE for EWMA for new cases with alpha=" + str(alpha)+ " is " + str(mse(error_hat_cases)))
  print("MAPE for EWMA for new deaths with alpha=" + str(alpha)+ " is " + str(mape(Y_obs_deaths[1:] , error_hat_deaths)) + "%")
  print("MSE for EWMA for new deaths with alpha=" + str(alpha)+ " is " + str(mse(error_hat_deaths)))


# + id="gwBg-9gWkELx" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 289} outputId="123825ad-fd60-486c-945f-78f95ac8268c" executionInfo={"status": "ok", "timestamp": 1589318300319, "user_tz": 240, "elapsed": 430, "user": {"displayName": "<NAME>", "photoUrl": "", "userId": "16347493411284187623"}}
def req_inf_1():
  fpath = 'us-states.csv'
  ar_3(current_df)
  ar_5(current_df)
  ewma(current_df,0.5)
  ewma(current_df,0.8)
req_inf_1()

# + [markdown] id="7LTi6cIQkKlv" colab_type="text"
# **Inference**
#
# As we can observe the MAPE and MSE reported are quite high. Reason for this is that the number of cases and deaths were growing somewhat exponentially for the dates under consideration.Since, AR(3) and AR(5) are regression models much suited for linear correlations, the **exponential curve** yield high error.
# EWMA  improves on simple variance by assigning weights to the periodic returns. By doing this, we can give greater weight to more recent returns. Hence error in EWMA smoothening is not very efficient due to the **exponential growth**. What happens here is that the weighted average is taken from previous true and estimated value i.e Error is being minimized given the previous data growth, however the new data point exceeds that error value and thus the smoothening doesn't happen easily. 
# We can also observe that the EWMA with $\alpha = 0.8 $ has less error as the weigthage for the true value is higher and the estimate for every new day is predicted closer to the previous day's actual value.

# + [markdown] id="CGdsLkSnX55n" colab_type="text"
# # Required Inferences 2 (Wald's, Z, T tests)

# + [markdown] id="6bGtOhu_iIVp" colab_type="text"
# Here onwards, when we mention null hypothesis, we  mean that last week's mean = second last week's mean. When we talk about mean, we are talking about number of new deaths and new cases per day. <br>
# In all 1 sample and 2 sample, previous week's mean would refer to different things in accordance with the questions. Basically in 1 sample tests, we take last week's mean as MLE of last week's mean. In 2 sample test, we take both the distribution's together and work on them as required in respective tests.

# + id="PY4Qi1hVL8g1" colab_type="code" outputId="db979f10-0f52-4cef-dae1-bc7517e6b8c6" colab={"base_uri": "https://localhost:8080/", "height": 289} executionInfo={"status": "ok", "timestamp": 1589314615734, "user_tz": 240, "elapsed": 531, "user": {"displayName": "<NAME>", "photoUrl": "", "userId": "16347493411284187623"}}
## DATA PREPARATION

curr_wk_end_dt = current_df.date.max()
curr_wk_st_dt = current_df.date.max()-timedelta(6)

prev_wk_st_dt = current_df.date.max() - timedelta(13)
prev_wk_end_dt = current_df.date.max() - timedelta(7)

# Last available week of data 
curr_filter = (current_df['date']>=curr_wk_st_dt) & (current_df['date']<=curr_wk_end_dt)
curr_wk_df = current_df.loc[curr_filter]

# Second to last available week of data
prev_filter = (current_df['date']>=prev_wk_st_dt) & (current_df['date']<=prev_wk_end_dt)
prev_wk_df = current_df.loc[prev_filter]

# alpha values
z_half_alpha = 1.96
z_alpha = 1.64
t_12_half_alpha = 2.1788
t_6_half_alpha = 2.4469

print(curr_wk_df)
print(prev_wk_df)


# + colab_type="code" outputId="83e595ea-5fdd-498a-c118-43b505ef651e" id="MjoMVEUXaCtQ" colab={"base_uri": "https://localhost:8080/", "height": 51} executionInfo={"status": "ok", "timestamp": 1589312535380, "user_tz": 240, "elapsed": 501, "user": {"displayName": "<NAME>", "photoUrl": "", "userId": "16347493411284187623"}}
# Wald's 1 population test 
def get_W_one_pop(prev_wk, curr_wk):
    theta_not = prev_wk.sum()/prev_wk.shape[0]
    theta_hat = curr_wk.sum()/curr_wk.shape[0]
    se_theta_hat = (np.sum(np.square(np.array(curr_wk) - theta_hat))/(curr_wk.shape[0]**2))**0.5
    w_stat = abs((theta_hat - theta_not)/se_theta_hat)
    ret_str = "H_0 rejected" if w_stat > z_half_alpha else "Failed to reject H_0"
    return w_stat, ret_str

print("W-statistic (mean of last 2 week's death count - 1 pop): "+ str(get_W_one_pop(prev_wk_df.deaths, curr_wk_df.deaths)))
print("W-statistic (mean of last 2 week's cases count - 1 pop): "+ str(get_W_one_pop(prev_wk_df.cases, curr_wk_df.cases)))


# + [markdown] id="oG6h_ie0aZ2S" colab_type="text"
# We can see that for number of cases in Wald's test, the hypothesis can't be rejected. In cases of deaths, we see that hypothesis is rejected. <br>
# Wald's test requirements aren't met i.e. we need AN estimator of mu. As our data sample is very small, we can't call it AN estimator of mu.

# + id="-Y58OuMCL8g3" colab_type="code" outputId="6cba4812-081b-4895-bf2d-0a9ee128c99a" colab={"base_uri": "https://localhost:8080/", "height": 50}
#Wald's 2 population test comparing means of last 2 weeks against each other
def get_W_two_pop(prev_wk, curr_wk):
    prev_wk_mean_mle = prev_wk.sum()/prev_wk.shape[0]
    curr_wk_mean_mle = curr_wk.sum()/curr_wk.shape[0]
    delta_hat = prev_wk_mean_mle - curr_wk_mean_mle
    # numerator of 2 population W statistic
    se_prev_component = np.sum(np.square(np.array(prev_wk) - prev_wk_mean_mle))/(prev_wk.shape[0]**2)
    se_curr_component = np.sum(np.square(np.array(curr_wk) - curr_wk_mean_mle))/(curr_wk.shape[0]**2)
    # denominator of 2 population W statistic
    se_hat = (se_prev_component + se_curr_component) ** 0.5
    w_statistic = abs(delta_hat/se_hat)
    ret_str = "H_0 rejected" if w_statistic > z_half_alpha else "Failed to reject H_0"
    return w_statistic, ret_str

print("W-statistic (mean of last 2 week's death count - 2 population): "+ str(get_W_two_pop(prev_wk_df.deaths, curr_wk_df.deaths)))
print("W-statistic (mean of last 2 week's cases count - 2 population): "+ str(get_W_two_pop(prev_wk_df.cases, curr_wk_df.cases)))


# + [markdown] id="lt2HWEu8eTjT" colab_type="text"
# We can see that for number of cases in Wald's 2 population test, the hypothesis can't be rejected. For number of deaths, we see that hypothesis is rejected. <br>
# Wald's test requirements aren't met i.e. we need AN estimator of mu. As both of our data samples (last and second last week numbers) are very small, we can't call MLE estimator as AN estimator of mu.

# + id="XTiyTpKwL8g9" colab_type="code" outputId="bb1da7c0-cb51-4983-cc75-1b4bb428a246" colab={"base_uri": "https://localhost:8080/", "height": 50}
# Z test 1 population
def get_zstat(prev_wk, curr_wk):
    theta_not = prev_wk.sum()/prev_wk.shape[0]
    sample_mean = curr_wk.sum()/curr_wk.shape[0]
    sample_sd = (np.sum(np.square(np.array(curr_wk) - sample_mean))/curr_wk.shape[0])**0.5
    z_stat = abs((sample_mean - theta_not)/(sample_sd/(curr_wk.shape[0]**0.5)))
    ret_str = "H_0 rejected" if z_stat > z_half_alpha else "Failed to reject H_0"
    return z_stat, ret_str

print("Z-statistic (mean of last 2 week's death count - 1 pop): "+ str(get_zstat(prev_wk_df.deaths, curr_wk_df.deaths)))
print("Z-statistic (mean of last 2 week's cases count - 1 pop): "+ str(get_zstat(prev_wk_df.cases, curr_wk_df.cases)))


# + [markdown] id="nZLVN7P_erS0" colab_type="text"
# For Z test, we see that hypothesis is rejected for number of death counts and we can't reject the hypothesis for number of cases. <br>
# We can't use z-test because of following reasons: <br>
# 1.   True standard deviation is not known
# 2.   n is small and we can't assume that data is normally distributed.
#
#

# + id="oySqCooxL8g_" colab_type="code" outputId="55e096bc-4a1d-4354-df8c-6969122092d3" colab={"base_uri": "https://localhost:8080/", "height": 50}
def get_tstat_one_pop(prev_wk, curr_wk):
    mu_not = prev_wk.sum()/prev_wk.shape[0]
    sample_mean = curr_wk.sum()/curr_wk.shape[0]
    sample_sd = ((np.sum(np.square(np.array(curr_wk) - sample_mean))/curr_wk.shape[0]))**0.5
    t_stat = abs((sample_mean - mu_not)/(sample_sd/(curr_wk.shape[0]**0.5)))
    ret_str = "H_0 rejected" if t_stat > t_6_half_alpha else "Failed to reject H_0"
    return t_stat, ret_str

print("T-statistic (mean of last 2 week's death count - 1 pop): "+ str(get_tstat_one_pop(prev_wk_df.deaths, curr_wk_df.deaths)))
print("T-statistic (mean of last 2 week's cases count - 1 pop): "+ str(get_tstat_one_pop(prev_wk_df.cases, curr_wk_df.cases)))


# + [markdown] id="Sk6T8scJfQvK" colab_type="text"
# For 1 population, we can't be sure if T-test is applicable because of following reasons: <br>
# 1.   We don't know if data is normally distributed
#
# However, we see that for death count, we see that hypothesis is rejected. For case count, we see that we fail to reject the hypothesis.
#

# + id="zNojSvInL8hB" colab_type="code" outputId="9fa5eecb-049e-4d16-b2bb-1220e8bfbd72" colab={"base_uri": "https://localhost:8080/", "height": 51} executionInfo={"status": "ok", "timestamp": 1589314172022, "user_tz": 240, "elapsed": 614, "user": {"displayName": "<NAME>", "photoUrl": "", "userId": "16347493411284187623"}}
def paired_ttest(prev_wk, curr_wk):
    d_dist = np.array(prev_wk) - np.array(curr_wk)
    d_mean = np.sum(d_dist)/d_dist.shape[0]
    d_sd = (np.sum(np.square(d_dist - d_mean))/d_dist.shape[0])**0.5
    t_stat = abs(d_mean/(d_sd/(d_dist.shape[0]**0.5)))
    ret_str = "H_0 rejected" if t_stat > t_6_half_alpha else "Failed to reject H_0"
    return t_stat, ret_str

print("Paired T-statistic (mean of last 2 week's death count - 2 pop): "+ str(paired_ttest(prev_wk_df.deaths, curr_wk_df.deaths)))
print("Paired T-statistic (mean of last 2 week's cases count - 2 pop): "+ str(paired_ttest(prev_wk_df.cases, curr_wk_df.cases)))


# + [markdown] id="eR9PGN-xfzeP" colab_type="text"
# For paired t-test, we see that null hypothesis is rejected in case of death counts and we can't reject null hypothesis for count of cases. <br>
# We don't know the distribution of D (i.e. difference of second last week's numbers and last week's numbers). We assume that it is normally distributed. Thus we can't be sure that paired ttest is applicable or not. <br>
# Degree of freedom is taken as 6 = 7-1.

# + id="2CRwsKqfL8hD" colab_type="code" outputId="b8eced58-7dbe-4b05-8f36-96f8ba2dd429" colab={"base_uri": "https://localhost:8080/", "height": 50}
def unpaired_ttest(prev_wk, curr_wk):
    prev_mean = np.sum(np.array(prev_wk))/prev_wk.shape[0]
    curr_mean = np.sum(np.array(curr_wk))/curr_wk.shape[0]
    d_mean = prev_mean - curr_mean
    var_prev = np.sum(np.square(np.array(prev_wk) - prev_mean))/prev_wk.shape[0]
    var_curr = np.sum(np.square(np.array(curr_wk) - curr_mean))/curr_wk.shape[0]
    t_stat = abs(d_mean/(((var_prev/prev_wk.shape[0])+(var_curr/curr_wk.shape[0]))**0.5))
    ret_str = "H_0 rejected" if t_stat > t_12_half_alpha else "Failed to reject H_0"
    return t_stat, ret_str
    
print("Unpaired T-statistic (mean of last 2 week's death count - 2 pop): "+ str(unpaired_ttest(prev_wk_df.deaths, curr_wk_df.deaths)))
print("Unpaired T-statistic (mean of last 2 week's cases count - 2 pop): "+ str(unpaired_ttest(prev_wk_df.cases, curr_wk_df.cases)))


# + [markdown] id="Ffxq_9yHgcIh" colab_type="text"
# Unpaired ttest assumes that both samples are normally distributed and they are independent. We know that the samples aren't independent, as number of new cases/deaths does indirectly impact next week's new cases and deaths. Thus, we can't apply unpaired ttest. <br>
# We see that for deaths, null hypothesis is rejected. For cases, we fail to reject the null hypothesis. <br>
# Degree of freedom is taken as 12 = 14-1-1.

# + [markdown] id="7YdAE5jjkZzM" colab_type="text"
# # Required Inference 3 - Permutation test and KS test

# + [markdown] id="kXMbPPoElKFb" colab_type="text"
# Hereonwards, when we mention null hypothesis, we mean that the two samples come from the same distribution.

# + id="snwzNKDbL8hF" colab_type="code" outputId="8f01e86a-8bc3-4855-fdf9-7fd307e13469" colab={"base_uri": "https://localhost:8080/", "height": 51} executionInfo={"status": "ok", "timestamp": 1589315465446, "user_tz": 240, "elapsed": 2035, "user": {"displayName": "<NAME>", "photoUrl": "", "userId": "16347493411284187623"}}
# Permutation test
def permutation_test(X, Y, n, threshold, metadata):
    x, y = sum(X)/len(X), sum(Y)/len(Y)
    len_x, len_y = len(X), len(Y)
    cnt = 0
    
    L = []
    L.extend(X)
    L.extend(Y)
    # We are taking random n permutations for finding pvalue
    for _ in range(n):
        Lp = np.random.permutation(L)
        xp, yp = sum(Lp[:len_x])/len_x, sum(Lp[len_x:])/len_y
        if abs(xp - yp) > abs(x - y):
            cnt += 1
    result = "H_0 rejected" if cnt/n<=threshold else "Failed to reject H_0"
    print("P-value ("+metadata+") with n="+str(n)+": "+str(cnt/n)+"  "+result)

permutation_test(prev_wk_df.deaths.tolist(), curr_wk_df.deaths.tolist(), 40000, 0.05, "deaths")
permutation_test(prev_wk_df.cases.tolist(), curr_wk_df.cases.tolist(), 40000, 0.05, "cases")


# + [markdown] id="1_gj9CSKhF78" colab_type="text"
# Permutation test is applicable as it has no assumptions. According to permutation test, we fail to reject that distribution of 2 samples are same. But in cases of deaths, we reject it. <br>
# We take 40000 permutations to test the hypothesis as taking all 14! permutations isn't computationally viable.

# + id="0IgK75orL8hI" colab_type="code" colab={}
# K-S test for 2 population
def KS_test_2_pop(X, Y, X_name, Y_name, threshold):
    X_list = sorted(X)
    Y_list = sorted(Y)
    if X_list[-1]>Y_list[-1]:
        X_list, Y_list = Y_list, X_list
        X_name, Y_name = Y_name, X_name
    X_cdf_left = np.array([i for i in range(0, len(X_list))])/len(X_list)
    X_cdf_right = np.array([i for i in range(1, len(X_list)+1)])/len(X_list)
    Y_cdf_left = np.array([i for i in range(0, len(Y_list))])/len(Y_list)
    Y_cdf_right = np.array([i for i in range(1, len(Y_list)+1)])/len(Y_list)
    
    # Table creation for final printing
    X_table = pd.DataFrame(X_list, columns=['data'])
    X_table['cdf_left'] = X_cdf_left
    X_table['cdf_right'] = X_cdf_right
    
    Y_table = pd.DataFrame(Y_list, columns=['data'])
    Y_table['cdf_left'] = Y_cdf_left
    Y_table['cdf_right'] = Y_cdf_right
    
    y_it = 0
    left_diff = []
    right_diff = []
    max_diff, x, y1, y2  = -math.inf, -1, -1, -1

    # We have many cases in the below for loop, covering all cases 
    # ranging from when both distribution have same data point and
    # when they have different data points. We also cover both sides
    # of a ecdf. left of ecdf and right of ecdf.
    for i in range(X_table.shape[0]):
        while Y_table.data[y_it]<=X_table.data[i]: y_it+=1
        if Y_table.data[y_it] == X_table.data[i]:
            left_diff.append(abs(X_table.cdf_left[i]-Y_table.cdf_left[y_it]))
            right_diff.append(abs(X_table.cdf_right[i]-Y_table.cdf_right[y_it]))
            if left_diff[-1]>max_diff:
                max_diff = left_diff[-1]
                x = X_table.data[i]
                y1 = min(X_table.cdf_left[i],Y_table.cdf_left[y_it])
                y2 = max(X_table.cdf_left[i],Y_table.cdf_left[y_it])
            elif right_diff[-1]>max_diff:
                max_diff = right_diff[-1]
                x = X_table.data[i]
                y1 = min(X_table.cdf_right[i],Y_table.cdf_right[y_it])
                y2 = max(X_table.cdf_right[i],Y_table.cdf_right[y_it])
        else:
            left_diff.append(abs(X_table.cdf_left[i]-Y_table.cdf_left[y_it]))
            right_diff.append(abs(X_table.cdf_right[i]-Y_table.cdf_left[y_it]))
            if left_diff[-1]>max_diff:
                max_diff = left_diff[-1]
                x = X_table.data[i]
                y1 = min(X_table.cdf_left[i],Y_table.cdf_left[y_it])
                y2 = max(X_table.cdf_left[i],Y_table.cdf_left[y_it])
            if right_diff[-1]>max_diff:
                max_diff = right_diff[-1]
                x = X_table.data[i]
                y1 = min(X_table.cdf_right[i],Y_table.cdf_left[y_it])
                y2 = max(X_table.cdf_right[i],Y_table.cdf_left[y_it])
    X_table['left_diff'] = left_diff
    X_table['right_diff'] = right_diff
    
    # calculations for drawing
    X_print = [[min(X_list[0], Y_list[0])-2, 0]]
    for i in range(len(X_list)):
        X_print.append([X_list[i], X_cdf_left[i]])
        X_print.append([X_list[i], X_cdf_right[i]])
    X_print.append([max(X_list[-1], Y_list[-1])+2, 1])
    X_print = pd.DataFrame(X_print)
    
    Y_print = [[min(X_list[0], Y_list[0])-2, 0]]
    for i in range(len(Y_list)):
        Y_print.append([Y_list[i], Y_cdf_left[i]])
        Y_print.append([Y_list[i], Y_cdf_right[i]])
    Y_print = pd.DataFrame(Y_print)
    
    fig, ax = plt.subplots(figsize=(8,6))
    plt.plot(X_print[0], X_print[1], color='red', label=X_name)
    plt.plot(Y_print[0], Y_print[1], color='blue', label=Y_name)
    plt.plot([x, x], [y1, y2], color='green', label='Max difference')
    ax.set_xlabel('Counts')
    ax.set_ylabel('eCDF')
    ax.set_title('K-S Test to find statistic (max cdf difference)')
#     ax.annotate("K-stat:" + str(round(max_diff,3)), xy=(x*1.005, min(y1,y2)+(y2-y1)/2))
    ax.text(x/max(max(X_list), max(Y_list))-0.1, min(y1,y2)+(y2-y1)/2, "K-stat: " + str(round(max_diff,3)), bbox={'facecolor': 'white', 'alpha': 0.8, 'pad': 10}, color='black', transform=ax.transAxes)
    ax.text(0.8, 0.1, 'H_0 rejected' if max_diff>threshold else 'H_1 accepted', bbox={'facecolor': 'white', 'alpha': 1, 'pad': 10}, color='red', transform=ax.transAxes)
    plt.legend()
    plt.show()
    print(X_name + " Table: ")
    print(X_table)



# + [markdown] id="3jNzvuABg8ye" colab_type="text"
# KS test is applicable as it has no assumptions.

# + id="4m2j6wrzL8hK" colab_type="code" outputId="94cdadec-4f0d-4ee3-ae35-d09ab7e18725" colab={"base_uri": "https://localhost:8080/", "height": 557} executionInfo={"status": "ok", "timestamp": 1589315617337, "user_tz": 240, "elapsed": 1096, "user": {"displayName": "<NAME>", "photoUrl": "", "userId": "16347493411284187623"}}
KS_test_2_pop(curr_wk_df.deaths.tolist(), prev_wk_df.deaths.tolist(), 'curr', 'prev', 0.05)


# + id="wjr7I7AfL8hN" colab_type="code" outputId="04548e58-353e-4fec-ddb9-b39486e748ff" colab={"base_uri": "https://localhost:8080/", "height": 557} executionInfo={"status": "ok", "timestamp": 1589315618901, "user_tz": 240, "elapsed": 642, "user": {"displayName": "<NAME>", "photoUrl": "", "userId": "16347493411284187623"}}
KS_test_2_pop(curr_wk_df.cases.tolist(), prev_wk_df.cases.tolist(), 'curr', 'prev', 0.05)


# + [markdown] id="p5Ds-1Fgl-xA" colab_type="text"
# The following are supporting functions for running 1 population K-S test.

# + id="jhRDNyQuL8hP" colab_type="code" colab={}
def factorial(number):
    if number == 0: return 1
    retVal = 1
    while number!=1:
        retVal*=number
        number-=1
    return retVal

# returns cdf of a binomail distribution at a point
def bin_cdf(param):
    x = param['x']
    n = param['n']
    p = param['p']
    return binom.cdf(x, n, p)
#     n_fact = factorial(n)
#     curr_tot = 0
#     for i in range(0, x+1):
#         curr_tot += ((1/(factorial(i)*factorial(n-i))) * (p**i) * ((1-p)**(n-i)))
#     return curr_tot * n_fact

# finds parameters of a binomial distribution given a sample
def get_bin_parameters(X):
    X_np = np.array(X)
    mean = np.sum(X_np)/X_np.shape[0]
    variance = np.sum(np.square(X_np - mean))/X_np.shape[0]
    p_hat = 1-variance/mean
    n_hat = mean**2 /(mean - variance)
    return n_hat, p_hat


# + id="-Bq1ZZ5XL8hR" colab_type="code" colab={}
# finds parameters of a geometric distribution given a sample
def get_geo_parameters(X):
    X_np = np.array(X)
    p_hat = X_np.shape[0]/np.sum(X_np)
    return p_hat

# returns cdf of a geometric distribution at a point
def geo_cdf(param):
    x = param['x']
    p = param['p']
    return geom.cdf(x, p)
#     curr_sum = 0
#     for i in range(1, x+1):
#         curr_sum += (1-p)**(i-1)
#     return p*curr_sum


# + id="Jj0yDm7HL8hT" colab_type="code" colab={}
# finds parameters of a poisson distribution given a sample
def get_poisson_parameters(X):
    X_np = np.array(X)
    lambda_hat = np.sum(X_np)/X.shape[0]
    return lambda_hat

# returns cdf of a poisson distribution at a point
def poisson_cdf(param):
    x = param['x']
    lmbda = param['lmbda']
    return poisson.cdf(x, lmbda)
#     curr_sum = 0
#     for i in range(0, x+1):
#         curr_sum += ((lmbda**i/math.factorial(i))*math.exp(-lmbda))
#     return curr_sum


# + id="_F2Ga7JXL8hV" colab_type="code" colab={}
# K-S test for 1 population
def KS_test_1_pop(X, function, parameters, X_name, threshold):
    X_list = sorted(X)
    X_cdf_left = np.array([i for i in range(0, len(X_list))])/len(X_list)
    X_cdf_right = np.array([i for i in range(1, len(X_list)+1)])/len(X_list)
    
    X_table = pd.DataFrame(X_list, columns=['data'])
    X_table['cdf_left'] = X_cdf_left
    X_table['cdf_right'] = X_cdf_right    

    left_diff = []
    right_diff = []
    max_diff, x, y1, y2  = -math.inf, -1, -1, -1
    for i in range(X_table.shape[0]):
        parameters['x'] = X_table.data[i]
        left_diff.append(abs(X_table.cdf_left[i] - function(parameters)))
        right_diff.append(abs(X_table.cdf_right[i] - function(parameters)))
        if left_diff[-1]>max_diff:
            max_diff = left_diff[-1]
            x = X_table.data[i]
            y1 = X_table.cdf_left[i]
            y2 = function(parameters)
        if right_diff[-1]>max_diff:
            max_diff = right_diff[-1]
            x = X_table.data[i]
            y1 = X_table.cdf_right[i]
            y2 = function(parameters)
    X_table['left_diff'] = left_diff
    X_table['right_diff'] = right_diff
    
    if y1>y2: y1,y2 = y2,y1

    # calculations for drawing
    distribution_print = []
    for i in range(min(X_table.data)-30, max(X_table.data)+31):
        parameters['x'] = i
        distribution_print.append([i, function(parameters)])
    distribution_print = np.array(distribution_print).transpose()
    
    fig, ax = plt.subplots(figsize=(8,6))
    plt.step([(min(X_table.data)-30)] + X_table.data.tolist() + [(max(X_table.data)+30)], [0]+ X_table.cdf_right.tolist()+[1], where='post', color='red', label=X_name)
    plt.plot(distribution_print[0], distribution_print[1], color='blue', label='Distribution')
    plt.plot([x, x], [y1, y2], color='green', label='Max difference')
    ax.set_xlabel('Counts')
    ax.set_ylabel('eCDF/CDF')
    ax.set_title('K-S Test to find statistic (max cdf difference)')
#     ax.annotate(round(max_diff,3), xy=(x*1.005, min(y1,y2)+(y2-y1)/2))
    ax.text(x/(max(X_table.data)+30)-0.1, min(y1,y2)+(y2-y1)/2, "K-stat: " + str(round(max_diff,3)), bbox={'facecolor': 'white', 'alpha': 0.8, 'pad': 10}, color='black', transform=ax.transAxes)
    ax.text(0.05, 0.75, 'H_0 rejected' if max_diff>threshold else 'H_1 accepted', bbox={'facecolor': 'white', 'alpha': 1, 'pad': 10}, color='red', transform=ax.transAxes)
    plt.legend()
    plt.show()
    print(X_name + " Table: ")
    print(X_table)



# + id="Ujjyk3bCL8hW" colab_type="code" outputId="bc05e17a-877f-47a6-9a93-19fe1267dfc2" colab={"base_uri": "https://localhost:8080/", "height": 555}
param = {'p': get_geo_parameters(prev_wk_df.deaths)} 
KS_test_1_pop(curr_wk_df.deaths.tolist(), geo_cdf, param, 'deaths_geometric', 0.05)

# + id="VmBB73aaL8hY" colab_type="code" outputId="43050462-2139-4b5e-89ab-b684dcaf628d" colab={"base_uri": "https://localhost:8080/", "height": 555}
param = {'p': get_geo_parameters(prev_wk_df.cases)} 
KS_test_1_pop(curr_wk_df.cases.tolist(), geo_cdf, param, 'cases_geometric', 0.05)

# + id="m_cPSTLcL8ha" colab_type="code" outputId="d412875d-7005-4bb4-d4e3-5c64a4f6d384" colab={"base_uri": "https://localhost:8080/", "height": 555}
param = {'lmbda': get_poisson_parameters(prev_wk_df.deaths)}
KS_test_1_pop(curr_wk_df.deaths.tolist(), poisson_cdf, param, 'deaths_poisson', 0.05)

# + id="ikKxu4iCL8hc" colab_type="code" outputId="477ffe11-2029-4d3b-e91e-791b93d6c4aa" colab={"base_uri": "https://localhost:8080/", "height": 555}
param = {'lmbda': get_poisson_parameters(prev_wk_df.cases)}
KS_test_1_pop(curr_wk_df.cases.tolist(), poisson_cdf, param, 'cases_poisson', 0.05)

# + id="PZinblbOL8he" colab_type="code" outputId="c686178a-877a-4cf5-a6a6-90de598bffd7" colab={"base_uri": "https://localhost:8080/", "height": 555}
n_hat, p_hat = get_bin_parameters(prev_wk_df.cases)
param = {'n': n_hat, 'p': p_hat}
KS_test_1_pop(curr_wk_df.cases.tolist(), bin_cdf, param, 'cases_binomial', 0.05)

# + id="txCHaKVjL8hg" colab_type="code" outputId="de564c77-26d7-418e-ef41-d5fb65dd7ada" colab={"base_uri": "https://localhost:8080/", "height": 555}
n_hat, p_hat = get_bin_parameters(prev_wk_df.deaths)
param = {'n': n_hat, 'p': p_hat}
KS_test_1_pop(curr_wk_df.deaths.tolist(), bin_cdf, param, 'deaths_binomial', 0.05)

# + [markdown] id="FhjhsSQUrH3S" colab_type="text"
# # Pearson Correlation - Required Inference 4

# + [markdown] id="7eH0wys7rPMD" colab_type="text"
# Here, we are going to observe the correlation of Closing price of IT Stocks with the #cases and #deaths of Covid Respectively.

# + id="oX5Cd8EnrK9h" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 34} outputId="805b194c-01b8-4b63-adc7-047e3fafa690" executionInfo={"status": "ok", "timestamp": 1589317004024, "user_tz": 240, "elapsed": 422, "user": {"displayName": "<NAME>", "photoUrl": "", "userId": "16347493411284187623"}}
import math  

correlation = []

Y_mean = df['Close'].mean()
Y_new_sum = ((df['Close'].apply(lambda x: x - Y_mean)).apply(lambda x: x*x)).sum()

X = ['cases', 'deaths']
for i in range(len(X)):
  X_mean = df[X[i]].mean()
  X_new_sum = ((df[X[i]].apply(lambda x: x - X_mean)).apply(lambda x: x*x)).sum()
  XY_sum = ((df[X[i]].apply(lambda x: x - X_mean)) * (df['Close'].apply(lambda x: x - Y_mean))).sum()

  correlation.append( XY_sum / (math.sqrt(Y_new_sum) * math.sqrt(X_new_sum)))
  
print (correlation)

# + [markdown] id="Okmij6PnrRmp" colab_type="text"
# As we can observe from the result, both #cases and #deaths have correlation less than 0.5 with the closing price of stocks.
# Thus, IT sector Stocks have little to NO correlation with the #cases and #deaths of Covid.

# + [markdown] id="MKBpsa5OS1y1" colab_type="text"
# # BAYESIAN INFERENCE - Required inference 5

# + id="3tcamxZ0S8xY" colab_type="code" outputId="ae1b7de3-284a-4468-c9ce-11ea2d79f2a6" executionInfo={"status": "ok", "timestamp": 1589315802608, "user_tz": 240, "elapsed": 487, "user": {"displayName": "<NAME>", "photoUrl": "", "userId": "16347493411284187623"}} colab={"base_uri": "https://localhost:8080/", "height": 221}
curr_filter = (current_df['date']>"2020-03-31") & (current_df['date']<"2020-05-01")
april_df = current_df.loc[curr_filter]

data = april_df.deaths

print(april_df.head())
print(data.head())

# + [markdown] id="05f3z-KgTBz8" colab_type="text"
# ### Function to plot gamma distribution 
#
# *   We observed that the posterior becomes a gamma distribution. (with parameters ( $(\sum X_i + 1$, $ n + \frac{1}{\bar{X}}$))
#
# *   So we plotted the data with the number of new deaths, considering data of 7 days, 14 days, 21 days and 28 days.
# *  Also, MAP for gamma$(\alpha, \beta)$ is $\frac{\alpha}{\beta}$

# + id="ImMUGlx3TDu1" colab_type="code" colab={}
from scipy.stats import gamma
import numpy as np
import matplotlib.pyplot as plt
def plot_gamma(a=1, b=1, label="0"):
  # fig, ax = plt.subplots(1, 1)
  x = np.linspace(gamma.ppf(0.01, a, scale=1/b),
                      gamma.ppf(0.99, a, scale=1/b), 1000)
  # ax.plot(x, gamma.pdf(x, a, scale=1/b), label=label)
  plt.plot(x, gamma.pdf(x, a, scale=1/b), label=label)
  plt.legend()


# + id="Pd9zxHl7TF5b" colab_type="code" colab={}
def plot_pdf(data):
  X = sum(data)
  l = len(data)
  alpha = X + 1
  x_hat_first_week = sum(data[:7])/7
  beta = (l + 1/x_hat_first_week)
  plot_gamma(alpha, beta, str(l) + " map: " + str(int(alpha/beta)) )


# + id="5eXhItoZTGxm" colab_type="code" outputId="2778a77e-e1fb-488a-ed5f-3ef40de88463" executionInfo={"status": "ok", "timestamp": 1589315858885, "user_tz": 240, "elapsed": 823, "user": {"displayName": "<NAME>", "photoUrl": "", "userId": "16347493411284187623"}} colab={"base_uri": "https://localhost:8080/", "height": 390}
plt.figure(figsize=(10,6))
plot_pdf(data[:7])
plot_pdf(data[:14])
plot_pdf(data[:21])
plot_pdf(data[:28])
plt.title("Posterior for Gamma")
plt.show()

# + [markdown] id="BZzwIOApTc88" colab_type="text"
# Results for Bayesian
#
#
# *   We observe that the posterior follows a gamma distribution and as the data increases, it converges close to the mean.
# *   The MAP value and the end of 4 weeks is 640.

# + [markdown] id="ul5_vUodRmB-" colab_type="text"
# # Chi Square Test - Sample Inference 1
# We are testing whether Covid has contributed to an unseen volatility in the market
#
# *   Hypothesis: COVID and market volatility are independent
# *   We can use chisquare as it doesn't have assumptions
# *   This hypothesis is particularly useful for long term investors as wrong entry and exit points in a volatile market can hinder profits and make us buy a stock for undesirable price
# *   This is a significant analysis for intra day trading as a volatile market is the best market for them
#
# For our analysis, a market is considered to be volatile if it's index changes more than 2% in a day.

# + id="KaAqz5W7L8hi" colab_type="code" outputId="472cdeb8-9397-480d-9854-ed96139b4cfa" executionInfo={"status": "ok", "timestamp": 1589310601742, "user_tz": 420, "elapsed": 508, "user": {"displayName": "<NAME>", "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14Ghn2jGwld41IRFNRJDz6HF_pjoyPR38aG5aLx9A=s64", "userId": "17762003251198703139"}} colab={"base_uri": "https://localhost:8080/", "height": 101}
stocks_df = all_stock_df.copy().rename(columns={'date':'Date Time'})

# Get the % change in the index for the day
stocks_df["Change_p"] = abs(round(all_stock_df['Change']/all_stock_df['Close'] * 100,2))
# stocks_df["Change_p"] = abs(round(data_stocks['Change'].diff() * 100,2))
# print(stocks_df.head())
stocks_df['Date Time'] = pd.to_datetime(stocks_df["Date Time"])

# stocks_filter = ((stocks_df['Date Time']>="2020-03-01") & (stocks_df['Date Time']<"2020-05-01")) or ((stocks_df['Date Time']>="2019-03-01") & (stocks_df['Date Time']<"2019-05-01"))
post_stocks_df = stocks_df[((stocks_df['Date Time']>="2020-03-01") & (stocks_df['Date Time']<"2020-05-01"))]
pre_stocks_df = stocks_df[((stocks_df['Date Time']>="2019-03-01") & (stocks_df['Date Time']<"2019-05-01"))]
stocks_df = stocks_df[((stocks_df['Date Time']>="2020-03-01") & (stocks_df['Date Time']<"2020-05-01")) | 
                      ((stocks_df['Date Time']>="2019-03-01") & (stocks_df['Date Time']<"2019-05-01"))]

def apply_map(x):
    if x>=pd.to_datetime("2020-03-01") and x<pd.to_datetime("2020-05-01"):
        return "post covid"
    else:
        return "pre covid"

def apply_map_volatility(x):
    if x>=2:
        return "volatile"
    else:
        return "not volatile"

stocks_df["Label"] = stocks_df["Date Time"].apply(apply_map)
stocks_df["Volatility"] = stocks_df["Change_p"].apply(apply_map_volatility)
stocks_df["Count"] = [1] * len(stocks_df["Change_p"])
# print(stocks_df.head())
stocks_df = stocks_df[["Label", "Volatility", "Count"]]

table = pd.pivot_table(stocks_df, values='Count', index=['Label'], columns=['Volatility'], aggfunc=np.count_nonzero)
print(table)
obs = [[table.loc["post covid", "not volatile"], table.loc["post covid", "volatile"]], 
     [table.loc["pre covid", "not volatile"], table.loc["pre covid", "volatile"]]]

g, p, dof, expctd = chi2_contingency(obs)
print("test statistic: ", g, "p value: ", p)
# df = stocks_df.set_index('Date Time').join(current_df.set_index('date'))
# print(df.head())

# + [markdown] id="6RkzLlCnSu0p" colab_type="text"
# Results for Sample Inference 1
#
# As you can see above, the p value is 10^-7 which is quiet low even for a 99% confidence.
#
# Hence we reject the hypothesis of Volatility and  covid being independent i.e, market volatility is highly dependent on covid 

# + [markdown] id="FEzJAgAhnmpH" colab_type="text"
# # Sample Inference 2
# As we see in below visualization that number of cases and deaths started increasing rapidly from March 1st 2020. Thus we wanted to analyze its impact on the stocks. <br>
# For that, we have taken samples for 35 days before March 1st (Sample precovid) and 35 days post March 1st (Sample postcovid). We then run Paired ttest on these two data samples. <br>
# We are running one-tailed test with null hypothesis as  mean of Sample precovid <= mean of sample postcovid. <br>

# + id="Nn7Vm1OcnKZW" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 391} outputId="32c128d0-5900-4a6b-a1b4-e3e991f66574" executionInfo={"status": "ok", "timestamp": 1589316515565, "user_tz": 240, "elapsed": 696, "user": {"displayName": "<NAME>", "photoUrl": "", "userId": "16347493411284187623"}}
fig, ax = plt.subplots(figsize=(15,6))
n=100
plt.plot(current_df.date.tolist()[:n], current_df.cases.tolist()[:n])
plt.plot(current_df.date.tolist()[:n], current_df.deaths.tolist()[:n])

# + id="Azc4n5tJL8hp" colab_type="code" colab={}
post_stocks_df = stock_df[stock_df['date']>="2020-03-01"]
pre_stocks_df = stock_df[stock_df['date']<"2020-03-01"]


# + id="Xi2ndFZwL8hz" colab_type="code" colab={}
def paired_ttest_one_tailed(pre_covid, post_covid):
    t_34_005 = 1.6909
    d_dist =  post_covid - pre_covid
    d_mean = np.sum(d_dist)/d_dist.shape[0]
    d_sd = (np.sum(np.square(d_dist - d_mean))/d_dist.shape[0])**0.5
    t_stat = d_mean/(d_sd/(d_dist.shape[0]**0.5))
    ret_str = "H_0 rejected" if t_stat < -t_34_005 else "Failed to reject H_0"
    return t_stat, ret_str



# + id="eEtcsIK-L8h1" colab_type="code" outputId="d47fd5a1-b657-49e0-8be9-eecfe737e5f8" colab={"base_uri": "https://localhost:8080/", "height": 34} executionInfo={"status": "ok", "timestamp": 1589316529751, "user_tz": 240, "elapsed": 443, "user": {"displayName": "<NAME>", "photoUrl": "", "userId": "16347493411284187623"}}
post_covid_stock_data = np.array(stock_df[stock_df['date']>="2020-03-01"].Close.tolist()[:35])
pre_covid_stock_data = np.array(stock_df[stock_df['date']<"2020-03-01"].Close.tolist()[-35:])

print("One tailed Paired T-statistic (mean of 35 days prior Covid and 35 days after Covid (03-01)): "+\
      str(paired_ttest_one_tailed(pre_covid_stock_data, post_covid_stock_data)))

# + [markdown] id="rN1Ck4cMpcRU" colab_type="text"
# We see that null hypothesis is rejected. **This means the increase in count of covid cases/deaths has negatively impacted the stocks market.**

# + [markdown] id="MeTO2OAarWuP" colab_type="text"
# # Sample Inference 3 - Lockdown Impacts

# + [markdown] id="48DxtfQpvXfo" colab_type="text"
# **Multi Linear regression**

# + id="tvDwkiC2sZLI" colab_type="code" colab={}
curr_filter = (covid_stock_df['date']>="2020-02-01") & (covid_stock_df['date']<"2020-03-20")
first_df = covid_stock_df.loc[curr_filter]
first_df
close_df=covid_stock_df.loc[curr_filter]

# + id="NaJFeFpPsayc" colab_type="code" colab={}
curr_filter1 = (covid_stock_df['date']>="2020-03-20") & (covid_stock_df['date']<="2020-04-30")
second_df = covid_stock_df.loc[curr_filter1]
close_second_df=covid_stock_df.loc[curr_filter1]

# + id="K56mC3qIsd0I" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 119} outputId="0a9325a4-2013-4074-b503-79629812ebde" executionInfo={"status": "ok", "timestamp": 1589317928399, "user_tz": 240, "elapsed": 908, "user": {"displayName": "<NAME>", "photoUrl": "", "userId": "16347493411284187623"}}
second_df['constant_value']=1 #Adding the constant  as 1

# + id="q0WbOlV8siU9" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 119} outputId="1b9381fb-2d63-4c5d-8a31-4d8fefe8a7a3" executionInfo={"status": "ok", "timestamp": 1589317928400, "user_tz": 240, "elapsed": 389, "user": {"displayName": "<NAME>", "photoUrl": "", "userId": "16347493411284187623"}}
first_df['constant_value']=1 #Adding the constant  as 1

# + id="mQ8ZV9QsslUS" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 374} outputId="1ded7037-89ce-4df2-eb09-88e306152f86" executionInfo={"status": "ok", "timestamp": 1589317931281, "user_tz": 240, "elapsed": 493, "user": {"displayName": "<NAME>", "photoUrl": "", "userId": "16347493411284187623"}}
y = np.array(first_df['Close'])
y_second=np.array(second_df['Close'])
del first_df['date']
del first_df['Change']
del second_df['date']
del second_df['Change']
cols = first_df.columns.tolist()
cols_second=second_df.columns.tolist()
cols = cols[-1:] + cols[:-1]
first_df=first_df[cols]
second_df=second_df[cols]
del first_df['Close']
cols = first_df.columns.tolist()
del second_df['Close']
cols = second_df.columns.tolist()
X=first_df.values
X_second=second_df.values
def fit(X,y):
  n = X.shape[0]
  P = X.shape[1] -1
  X_bar = np.mean(X,0)
  y_bar = np.mean(y)
  XX = X.T @ X
  beta = None
  y_hat= None
  SS_tot, SS_exp = None, None
  error = None
  beta = np.linalg.inv(XX) @ X.T @ y  #Calculating the coefficients
  y_hat = X @ beta #predicting the values
  prediction = y_hat
  residuals= y - y_hat # difference between predicted and actual values
  model=(X @ beta) + residuals # Multi Linear Regression model


  S = np.sum( (y - y_hat)**2 )
  SS_exp = np.sum( (y_hat - y_bar)**2 )
  SS_tot = np.sum( (y - y_bar)**2 )
  rsquarederror = SS_exp /SS_tot
  meansquareerror=S / (n - P -1)
  #S2_n_p_1 = S / (n - P -1)

  error=np.sqrt(meansquareerror)
  estimatesAndErrors=[]
  estimatesAndErrors.append(beta)
  estimatesAndErrors.append(prediction)
  estimatesAndErrors.append(residuals)
  output={'beta_values':beta,'predicted values':prediction,'residuals':residuals}
  print('-----ERROR-----------',error)
  return output
print("Linear regression for pre lock down scenario")
linearRegression=fit(X,y)
linearRegression

# + id="SYeW07s6uOPu" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 306} outputId="a56480ec-8f3e-4f1b-e719-06b6986d1bf9" executionInfo={"status": "ok", "timestamp": 1589317943888, "user_tz": 240, "elapsed": 966, "user": {"displayName": "<NAME>", "photoUrl": "", "userId": "16347493411284187623"}}
print("Linear regression for post lock down scenario")
linearRegrssion_second=fit(X_second,y_second) #applying linear regression for post lock down start data
linearRegrssion_second

# + id="W7P6fpL7uPyS" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 248} outputId="4f4eeb64-9647-4094-8875-c6fe1b9e0dca" executionInfo={"status": "ok", "timestamp": 1589317945895, "user_tz": 240, "elapsed": 1099, "user": {"displayName": "<NAME>", "photoUrl": "", "userId": "16347493411284187623"}}
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(first_df['cases'],first_df['deaths'],close_df['Close'],c='blue', marker='o', alpha=0.5)
ax.set_xlabel('cases')
ax.set_ylabel('deaths')
ax.set_zlabel('close')
plt.show()

# + id="7eA6pVuruVFV" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 248} outputId="40d40688-8fea-4f1e-f9b4-f5d7fd93074d" executionInfo={"status": "ok", "timestamp": 1589317947065, "user_tz": 240, "elapsed": 635, "user": {"displayName": "<NAME>", "photoUrl": "", "userId": "16347493411284187623"}}
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(second_df['cases'],second_df['deaths'],close_second_df['Close'],c='blue', marker='o', alpha=0.5)
ax.set_xlabel('cases')
ax.set_ylabel('deaths')
ax.set_zlabel('close')
plt.show()

# + [markdown] id="11qbVhlNuZJj" colab_type="text"
# On assumption that the relation between stock prices, number of deaths, number of cases is linear,we tried to estimate the coefficients of Multi Linear Regression model and predicted the stock prices pre lock down , post lock down. The data is divided in to two parts, one is from Feb 1 to March 19 and the other is  from March 20 to April 30. The lockdown in California started on March 20, Hence the division in the data to see the impact on Tech stocks. We calculated the error values and found that the RMSE is around 108 and 73. But from the 3D scatter plots between stocks,deaths,cases we can notice that relation is not linear which explains the reason for  high RMSE.Both during the pre-lock down and post lock down,Although there is some dependency, the stock prices aren't linearly dependent on no.of deaths and no.active cases based on 2 months data we analyzed. 

# + id="5D_QTziyudip" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 159} outputId="aec5686f-c818-41a1-9f9f-7f654a1c707e" executionInfo={"status": "ok", "timestamp": 1589317949666, "user_tz": 240, "elapsed": 982, "user": {"displayName": "<NAME>", "photoUrl": "", "userId": "16347493411284187623"}}
print('Parameter Estimates for data pre lock down')
    
df3 = pd.DataFrame(data = {'Parameter Estimate':linearRegression['beta_values']}, index = ['intercept', 'cases', 'deaths'])
        
display(df3.round(3))

# + id="yys8DbroueZB" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 159} outputId="8ad01d8d-f583-4497-e5b7-afbd7ef7bbed" executionInfo={"status": "ok", "timestamp": 1589317951208, "user_tz": 240, "elapsed": 540, "user": {"displayName": "<NAME>", "photoUrl": "", "userId": "16347493411284187623"}}
print('Parameter Estimates for data post lock down')
    
df3 = pd.DataFrame(data = {'Parameter Estimate':linearRegrssion_second['beta_values']}, index = ['intercept', 'cases', 'deaths'])
        
display(df3.round(3))

# + [markdown] id="JZDs_Y0lu9Jy" colab_type="text"
# **Pearson Correlation between pre-covid and post-covid stocks data**
# <br>
#

# + id="740eaVOAL8h7" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 34} outputId="4e9b0e78-25a8-4647-df93-844a2e797d37" executionInfo={"status": "ok", "timestamp": 1589317086279, "user_tz": 240, "elapsed": 826, "user": {"displayName": "<NAME>", "photoUrl": "", "userId": "16347493411284187623"}}
precovid = (all_stock_df['date']>="2019-02-01") & (all_stock_df['date']<="2019-05-01")
postcovid = (all_stock_df['date']>="2020-02-01") & (all_stock_df['date']<"2020-05-01")
stocks_precovid = all_stock_df.loc[precovid]
stocks_postcovid = all_stock_df.loc[postcovid]
# stocks_df = pd.concat([stocks_precovid, stocks_postcovid], axis=1, sort=False)
print (stocks_precovid.shape, stocks_postcovid.shape)

# + [markdown] id="J-aSdpQyvC1D" colab_type="text"
# This compares the correlation between stocks for the months of february-april for 2019 vs 2020 i.e. pre-covid Stocks' Closing Price vs Post-covid Stocks' Closing Price.

# + id="9xYOUqxArh_V" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 34} outputId="d7154bde-4618-4b6d-8532-4eab4f41202c" executionInfo={"status": "ok", "timestamp": 1589317086556, "user_tz": 240, "elapsed": 911, "user": {"displayName": "<NAME>", "photoUrl": "", "userId": "16347493411284187623"}}
import math  
import statistics 
X = stocks_precovid['Close'].values.tolist()
Y = stocks_postcovid['Close'].values.tolist()

Y_mean = statistics.mean(Y) 
Y_new = [number - Y_mean for number in Y]
Y_new_2 = [number*number for number in Y_new]
Y_new_sum = sum(Y_new_2)

X_mean = statistics.mean(X) 
X_new = [number - X_mean for number in X]
X_new_2 = [number*number for number in X_new]
X_new_sum = sum(X_new_2)

X_new_Y_new = [a*b for a, b in zip(X_new, Y_new)]
XY_sum = sum(X_new_Y_new)

correlation = XY_sum / (math.sqrt(Y_new_sum) * math.sqrt(X_new_sum))
  
print(correlation)

# + [markdown] id="bILYI_SavGPE" colab_type="text"
# From the above value, we know that the pre-covid and post-covid Stocks' Closing Price is negatively correlated as expected.

# + [markdown] id="gv-fqz7Zrp5T" colab_type="text"
# **Pearson Correlation between two date ranges of Covid**<br>
#
#
# 1.   Here, we have taken the first range from february 2020 to March 19, 2020. This signifies the period when Covid began but Lockdown had NOT been enforced in the any of the 3 states.
# 2.   The second range is from March 20, 2020 till April 30, 2020. Lockdown in California began on March 20, in New York shortly after that and in Washington just a few days before that.
#
#

# + id="6iRUQEkJrvBE" colab_type="code" colab={}
curr_filter = (covid_stock_df['date']>="2020-02-01") & (covid_stock_df['date']<"2020-03-20")
covid_start = covid_stock_df.loc[curr_filter]
curr_filter = (covid_stock_df['date']>="2020-03-20") & (covid_stock_df['date']<"2020-05-01")
covid_end = covid_stock_df.loc[curr_filter]

# + [markdown] id="Sr0MAh5trwGn" colab_type="text"
#
#
# 1.   From the first correlation between #cases and IT Stocks Closing Price we discover a very high negative correlation. This means as #cases increased, the closing price went down.
# 2.   The same is observed for correlation between #deaths and IT Stocks Closing Price. 
# 3.   This is the period when Lockdown had not been enforced.
#
#
#
#

# + id="lf-Z1T25ryPD" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 34} outputId="9afa736e-2c83-4b94-c3b5-79288b3ffff0" executionInfo={"status": "ok", "timestamp": 1589317158098, "user_tz": 240, "elapsed": 924, "user": {"displayName": "Bishal <NAME>", "photoUrl": "", "userId": "16347493411284187623"}}
import math  

correlation = []

Y_mean = covid_start['Close'].mean()
Y_new_sum = ((covid_start['Close'].apply(lambda x: x - Y_mean)).apply(lambda x: x*x)).sum()

X = ['cases', 'deaths']
for i in range(len(X)):
  X_mean = covid_start[X[i]].mean()
  X_new_sum = ((covid_start[X[i]].apply(lambda x: x - X_mean)).apply(lambda x: x*x)).sum()
  XY_sum = ((covid_start[X[i]].apply(lambda x: x - X_mean)) * (covid_start['Close'].apply(lambda x: x - Y_mean))).sum()

  correlation.append( XY_sum / (math.sqrt(Y_new_sum) * math.sqrt(X_new_sum)))
  
print (correlation)

# + [markdown] id="hXZtKmYlr0Wq" colab_type="text"
# 1.   From the first correlation between #cases and IT Stocks Closing Price we discover almost no correlation. This means the #cases had no effect on the closing price.
# 2.  For correlation between #deaths and IT Stocks Closing Price, we see a positive correlation. 
# 3.   This is the period when Lockdown had been enforced.
#

# + id="zFp6BN-br2lG" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 34} outputId="3dc60b7f-49a9-4e3b-8dab-b607829c7c28" executionInfo={"status": "ok", "timestamp": 1589317176318, "user_tz": 240, "elapsed": 918, "user": {"displayName": "<NAME>", "photoUrl": "", "userId": "16347493411284187623"}}
import math  

correlation = []

Y_mean = covid_end['Close'].mean()
Y_new_sum = ((covid_end['Close'].apply(lambda x: x - Y_mean)).apply(lambda x: x*x)).sum()

X = ['cases', 'deaths']
for i in range(len(X)):
  X_mean = covid_end[X[i]].mean()
  X_new_sum = ((covid_end[X[i]].apply(lambda x: x - X_mean)).apply(lambda x: x*x)).sum()
  XY_sum = ((covid_end[X[i]].apply(lambda x: x - X_mean)) * (covid_end['Close'].apply(lambda x: x - Y_mean))).sum()

  correlation.append( XY_sum / (math.sqrt(Y_new_sum) * math.sqrt(X_new_sum)))
  
print (correlation)

# + [markdown] id="Il6_OKBNr41Z" colab_type="text"
# 1. From the inferences observed above, we see that Lockdown, especially in California- the IT hub, had a great impact on the IT market. 
# 2. The pearson correlation in the required inferences also makes sense now, that for the overall period, there is no correlation.
# 3. This may be due to the sharp increase in negative correlation in the first half (before Lockdown) followed by increase in positive correlation in the second half (After lockdown) which balances out the overall effect to be no correlation.

# + id="N1p1yx-Fril5" colab_type="code" colab={}


StarcoderdataJupyterScriptsDedupFiltered