import numpy as np, pandas as pd
from IPython.display import display, Image, Latex, HTML
import statsmodels.api as sm; from statsmodels.iolib.summary2 import summary_col
import imageio, os

import informationalTable
import countryLevelGraphs
import cantonLevelGraphs
import residGraph

Welcome!

Created 4-24-2022 by Lucas Voichick

The following code shows the data processing pipeline that I created for my Honors Thesis.

You can click here to read my thesis

.

Next line runs only if data has not yet been concatted. Gzip version is imported in next cell

# Imports and concats all NEUZU data

if os.path.exists('NewRegistrationData.gz') == False: # If pre-concatted data is available, concat and gzip to a single dataframe
    assert False, "Wait! Comment this out only if you have the original paid data that includes PLZ information along with computer memory and time to spare. Otherwise just import the zipped csv file on the next line"
    Dataframe_data2016 = pd.read_excel('data/NEUZU-2016.xlsx', engine='openpyxl')
    Dataframe_data2017 = pd.read_excel('data/NEUZU-2017.xlsx', engine='openpyxl')
    Dataframe_data2018 = pd.read_excel('data/NEUZU-2018.xlsx', engine='openpyxl')
    Dataframe_data2019 = pd.read_excel('data/NEUZU-2019.xlsx', engine='openpyxl')
    Dataframe_data2020 = pd.read_excel('data/NEUZU-2020 delfirstline.xlsx', engine='openpyxl')
    Dataframe_data2021 = pd.read_excel('data/NEUZU-2021 delfirstline.xlsx', engine='openpyxl')

    #combining all data into a dataframe and saving as a .csv file and gzipped .csv file
    allNEUZU_together = pd.concat([Dataframe_data2016,Dataframe_data2017,Dataframe_data2018, Dataframe_data2019,Dataframe_data2020,Dataframe_data2021])
    allNEUZU_together = allNEUZU_together.drop(columns=['PLZ'])
    allNEUZU_together.to_csv('NewRegistrationData.gz', compression='gzip')

Import all data

allNEUZU_together = pd.read_csv('NewRegistrationData.gz', low_memory = False)
print("Imported entire NEUZU dataset, len:", "{:,}".format(len(allNEUZU_together)))

Imported entire NEUZU dataset, len: 2,425,320

Filtering, Aggregating at Quarter Level, Classifying by Electric / Non-Electric

### # # # # # # # # # # # # # # # # ##
## Filtering Unwanted Registrations ##
### # # # # # # # # # # # # # # # # ##

## removing car class =! 1, 2, 10 to keep only personal vehicles ##
carClassFilter = allNEUZU_together['Fahrzeugart_Code'].isin([1,2,10])
filteredNEUZU = allNEUZU_together[carClassFilter] 

## removing some anomalies with unexpected number of seats ##
filteredNEUZU = filteredNEUZU[~((filteredNEUZU['Sitzplätze'] <= 1.0) | (filteredNEUZU['Sitzplätze'] > 39.0))]

## removing those registered in unknown cantons ##
filteredNEUZU = filteredNEUZU[filteredNEUZU['Staat'] == 'Schweiz']

## removing registrations with missing datapoints ##
filteredNEUZU = filteredNEUZU.dropna(subset=['Erstinverkehrsetzung_Monat','Inverkehrsetzung_Kanton','Erstinverkehrsetzung_Jahr','Treibstoff', 'Treibstoff_Code'])

## removing datapoints with mislabelled canton registration ##
cantons_abbr = ['AG', 'AR', 'AI', 'BL', 'BS', 'BE', 'FR', 'GE', 'GL', 'GR', 'JU', 'LU', 'NE', 'NW', 'OW', 'SH', 'SZ', 'SO', 'SG', 'TG', 'TI', 'UR', 'VS', 'VD', 'ZG', 'ZH']
filteredNEUZU = filteredNEUZU[filteredNEUZU['Inverkehrsetzung_Kanton'].isin(cantons_abbr)]
print("Removed", "{:,}".format(len(allNEUZU_together) - len(filteredNEUZU)), "unwanted data points.", "{:,}".format(len(filteredNEUZU)), "remain")

Removed 692,600 unwanted data points. 1,732,720 remain

# Appending a quarter term
monthConditions = [
    (filteredNEUZU['Erstinverkehrsetzung_Monat'] <= 3),
    (filteredNEUZU['Erstinverkehrsetzung_Monat'] >= 4) & (filteredNEUZU['Erstinverkehrsetzung_Monat'] <= 6),
    (filteredNEUZU['Erstinverkehrsetzung_Monat'] >= 7) & (filteredNEUZU['Erstinverkehrsetzung_Monat'] <= 9),
    (filteredNEUZU['Erstinverkehrsetzung_Monat'] >= 10)]

filteredNEUZU['Quarter'] = np.select(monthConditions, [1,2,3,4])
filteredNEUZU['YearQuarter'] = filteredNEUZU['Erstinverkehrsetzung_Jahr'].astype(str)+'-'+filteredNEUZU['Quarter'].astype(str)
filteredNEUZU['YearQuarterCanton'] = filteredNEUZU['YearQuarter'] + '-' + filteredNEUZU['Inverkehrsetzung_Kanton']

# Separating by electric and non-electric
electricFilter = (filteredNEUZU['Treibstoff_Code'] == 'E') | (filteredNEUZU['Treibstoff_Code'] == 'X') | ((filteredNEUZU['Treibstoff_Code'] == 'R') & (filteredNEUZU['Hubraum'] == 0)) 
NEUZU_electric = filteredNEUZU[electricFilter]
NEUZU_nonEV    = filteredNEUZU[~electricFilter]

print("Of the", "{:,}".format(len(NEUZU_nonEV) + len(NEUZU_electric)), "registrations between 2016-2021,","{:,}".format(len(NEUZU_nonEV)),
      "are not fully electric vehicles and", "{:,}".format(len(NEUZU_electric)),"are fully electric")

Of the 1,732,720 registrations between 2016-2021, 1,654,831 are not fully electric vehicles and 77,889 are fully electric

Organizing data by canton, registration quarter-year, eletric vehicle, EV Ratio

Horizontal and Vertical Organization

# Breaks up electric and non-electric registrations by electric vehicle, non electric vehicle, 
# all registrations, and EV registration fraction at the canton level and returns them as a dataframe
def EV_NonEv_AllRegistrations_EVFration(canton):

    NotEVsNEUZU_byQuarterCanton = NEUZU_nonEV[NEUZU_nonEV['Inverkehrsetzung_Kanton'].isin(canton)]['YearQuarter'].value_counts(
        ).reindex(NEUZU_nonEV.YearQuarter.unique(), fill_value=0).sort_index(ascending=True).to_frame()
    
    EVsNEUZU_byQuarterCanton = NEUZU_electric[NEUZU_electric['Inverkehrsetzung_Kanton'].isin(canton)]['YearQuarter'].value_counts(
        ).reindex(NEUZU_electric.YearQuarter.unique(), fill_value=0).sort_index(ascending=True).to_frame()

    ## If all cantons are passed through as an arg, set name to 'All'
    canton = len(canton) > 1 and '_All' or "_"+ str(canton).replace("[","").replace("]","").replace("'","")
    
    YQRegistrations = pd.concat([EVsNEUZU_byQuarterCanton,NotEVsNEUZU_byQuarterCanton], axis = 1)
    YQRegistrations.columns = ['EVs' + canton, 'NonEVs' + canton]
    
    #calculating total number of vehicles, electric vehicle ratio
    YQRegistrations['AllVehicles'+canton] = YQRegistrations['EVs' + canton] + YQRegistrations['NonEVs'+canton]
    YQRegistrations['EVRatio'+canton] = (YQRegistrations['EVs' + canton] / ( YQRegistrations['AllVehicles'+canton] ) *100).round(4)
    return YQRegistrations

Horizontal Organization

#Calculates EV, nonEV, all, and Ratio for each canton, individually

YQRegistrations = pd.DataFrame([])
for canton in cantons_abbr:     
    YQRegistrations = pd.concat([YQRegistrations,EV_NonEv_AllRegistrations_EVFration([canton])], axis = 1) 
    
#calculates EV, nonEV, all, and Ratio for all cantons together
YQRegistrations = pd.concat([YQRegistrations,EV_NonEv_AllRegistrations_EVFration(cantons_abbr)], axis = 1).reset_index().rename(columns={"index": "YearQuarter"})

YQRegistrations.to_csv("EVData_HorizontallyOrganized.csv")

Horizontally organized df has all four Data Terms for each canton across all quarters

Sampled selection of row and columns shown below:

dataTerm   = 'EVs'#options include EVs, NonEVs, AllVehicles, EVRatio (see table at right)
cantonAbbr = 'BL' # 26 different canton options for graph shown in center table below (see table in center)

informationalTable.HorizontalData_InformationalGraphAndTable(dataTerm, cantonAbbr)
pd.options.display.max_columns = 100
YQRegistrations.iloc[1::2, 81:]

	EVs_TI	NonEVs_TI	AllVehicles_TI	EVRatio_TI	EVs_UR	NonEVs_UR	AllVehicles_UR	EVRatio_UR	EVs_VS	NonEVs_VS	AllVehicles_VS	EVRatio_VS	EVs_VD	NonEVs_VD	AllVehicles_VD	EVRatio_VD	EVs_ZG	NonEVs_ZG	AllVehicles_ZG	EVRatio_ZG	EVs_ZH	NonEVs_ZH	AllVehicles_ZH	EVRatio_ZH	EVs_All	NonEVs_All	AllVehicles_All	EVRatio_All
1	29	5030	5059	0.5732	1	321	322	0.3106	26	3661	3687	0.7052	59	8940	8999	0.6556	24	1697	1721	1.3945	216	13525	13741	1.5719	803	85761	86564	0.9276
3	33	5582	5615	0.5877	1	373	374	0.2674	26	3730	3756	0.6922	63	9269	9332	0.6751	27	1539	1566	1.7241	221	13150	13371	1.6528	955	83733	84688	1.1277
5	42	5207	5249	0.8002	4	380	384	1.0417	36	3555	3591	1.0025	92	9148	9240	0.9957	36	1862	1898	1.8967	210	13842	14052	1.4944	1008	85986	86994	1.1587
7	54	5414	5468	0.9876	4	302	306	1.3072	63	3526	3589	1.7554	113	9000	9113	1.2400	35	1744	1779	1.9674	505	13050	13555	3.7256	1627	80620	82247	1.9782
9	41	4917	4958	0.8269	1	341	342	0.2924	33	3173	3206	1.0293	95	8683	8778	1.0823	32	1741	1773	1.8049	288	13205	13493	2.1344	1165	85198	86363	1.3490
11	55	4621	4676	1.1762	1	222	223	0.4484	43	2527	2570	1.6732	90	7828	7918	1.1367	47	1478	1525	3.0820	433	11772	12205	3.5477	1700	72263	73963	2.2984
13	89	4834	4923	1.8078	4	314	318	1.2579	95	3137	3232	2.9394	319	8687	9006	3.5421	80	1650	1730	4.6243	650	12378	13028	4.9893	2870	83548	86418	3.3211
15	187	5112	5299	3.5290	8	242	250	3.2000	111	2991	3102	3.5783	518	8497	9015	5.7460	87	1559	1646	5.2855	971	12941	13912	6.9796	4266	80364	84630	5.0408
17	112	2886	2998	3.7358	8	199	207	3.8647	70	1839	1909	3.6668	312	4893	5205	5.9942	56	1027	1083	5.1708	556	7212	7768	7.1576	2531	46655	49186	5.1458
19	375	4407	4782	7.8419	26	235	261	9.9617	340	2490	2830	12.0141	763	6971	7734	9.8655	217	1412	1629	13.3211	2075	10389	12464	16.6479	8302	65458	73760	11.2554
21	400	3697	4097	9.7632	26	252	278	9.3525	364	2388	2752	13.2267	658	6657	7315	8.9952	186	1380	1566	11.8774	1383	9706	11089	12.4718	7516	62515	70031	10.7324
23	510	3046	3556	14.3420	43	181	224	19.1964	511	1942	2453	20.8316	1188	5147	6335	18.7530	283	1176	1459	19.3968	2305	7272	9577	24.0681	11288	47837	59125	19.0918

Organizes vertically for regression

Regression_CYQ = pd.DataFrame([])
mergedDf = pd.DataFrame([])

for canton in cantons_abbr:
    CantonVals = EV_NonEv_AllRegistrations_EVFration([canton])
    CantonVals['CantonName'] = canton
    temp = pd.concat([pd.DataFrame([]),CantonVals], axis = 0).reset_index().rename(columns={"index": "YearQuarter", "EVs_" + canton: "EVs","NonEVs_" + canton: "NonEVs","AllVehicles_" + canton: "AllVehicles","EVRatio_" + canton: "EVRatio"})
    Regression_CYQ = pd.concat([Regression_CYQ, temp], ignore_index = True)

Regression_CYQ['CantonYearQuarter'] = Regression_CYQ['CantonName'] + Regression_CYQ['YearQuarter'].astype(str)
Regression_CYQ['CantonYear'] = Regression_CYQ['CantonName'] +"-"+ Regression_CYQ['YearQuarter'].str[0:4]
Regression_CYQ['Year'] =    Regression_CYQ['YearQuarter'].str[0:4].astype(int)
Regression_CYQ['Quarter'] = Regression_CYQ['YearQuarter'].str[5].astype(int)

Appending regressors (purchase premiums, charging subsidies, tax ratio, population density)

# Appending a purchase premium term
purchasePremium = [
    (Regression_CYQ['CantonName'] == 'TG') & (Regression_CYQ['Year'] >= 2021),
    (Regression_CYQ['CantonName'] == 'TG') & (Regression_CYQ['Year'] == 2020),
    (Regression_CYQ['CantonName'] == 'TG') & (Regression_CYQ['Year'] == 2019),
    
    (Regression_CYQ['CantonName'] == 'TI') & ((Regression_CYQ['Year'] >= 2020) | ((Regression_CYQ['Year'] >= 2019) & (Regression_CYQ['Quarter'] >= 3))),
    
    (Regression_CYQ['CantonName'] == 'VS') & ((Regression_CYQ['Year'] == 2021) | ((Regression_CYQ['Year'] >= 2020) & (Regression_CYQ['Quarter'] == 4))),

    (Regression_CYQ['CantonName'] == 'SH') & ((Regression_CYQ['Year'] >= 2021))
]

premium_1000 = [2.0,3.5,4.0,  2,   3.5,   2]

Regression_CYQ['PP_1000'] = np.select(purchasePremium, premium_1000)
Regression_CYQ['PP_1000'] = np.where(Regression_CYQ['PP_1000'] > 0, Regression_CYQ['PP_1000'], 0)

Info on Purchase Premium Policies:

Thurgau
In 2019, a premium of 4000 CHF is offered
In 2020, decrease offer to 3500 CHF
In 2021, decrease to 2000 CHF

Ticino
Starting July 1, 2019, a premium of 2000 CHF is offered

Valais
From Nov. 1 2020 to end of 2022, a premium of 3500 CHF is offered

Schaffhausen
Starting in Jan 2021, a premium of 2000 CHF is offered

# Appending a charging station subsidy term
csSubsidy = [
    (Regression_CYQ['CantonName'] == 'GE') &  (Regression_CYQ['Year'] >= 2019),
    
    (Regression_CYQ['CantonName'] == 'TI') & ((Regression_CYQ['Year'] >= 2020) | ((Regression_CYQ['Year'] >= 2019) & (Regression_CYQ['Quarter'] >= 3))),
    
    (Regression_CYQ['CantonName'] == 'VS') & ((Regression_CYQ['Year'] == 2021) | ((Regression_CYQ['Year'] >= 2020) & (Regression_CYQ['Quarter'] == 4))),

    (Regression_CYQ['CantonName'] == 'VD') &  (Regression_CYQ['Year'] >= 2021),

]

cs_1000 = [1.0, 0.5, 1.5, 2.0]

Regression_CYQ['CS_1000'] = np.select(csSubsidy, cs_1000)
Regression_CYQ['CS_1000'] = np.where(Regression_CYQ['CS_1000'] > 0, Regression_CYQ['CS_1000'], 0) 

# Uncomment below to see canton-year-quarters that include either a purchase premium subsidy or charging station subsidy
# Regression_CYQ[(Regression_CYQ['PP_1000'] != 0.0) | (Regression_CYQ['CS_1000'] != 0.0)][['PP_1000', 'CS_1000', 'YearQuarter','CantonName']] 

Info on Charging Station Subsidies:

Geneva
Starting Jan 1 2029, a subsidy of 2000 CHF is offered

Ticino
Starting June 19, 2019, a subsidy of 500 CHF is offered (which was no longer offered by end of 2021)

Valais
Starting Nov 1 2020, a subsidy of 1500 CHF is offered (for a single 11-22 kW charging station)

Vaud
Starting Dec 2020, a subsidy of 2000 CHF is offered

(More info on charging station subsidies)

Appending Vehicle Price Ratio and Social Controls (Tertiary and Pop Density)

CPDY = pd.read_csv('Population&TertiaryData.gz')
CYQ_preDummies = pd.merge(left=Regression_CYQ, right=CPDY, left_on='CantonYear', right_on='Canton-Year').rename(columns={"Year_x": "Year"})

CYQ_preDummies['Pct Tertiary']       = CYQ_preDummies['Pct_Tertiary'] * 100
CYQ_preDummies['PopDen(Over25)_log'] = np.log(CYQ_preDummies['Population_Density_Over25'])

Info on controls:

Vehicle Price Ratio Controls: Touring Club Suisse - 2021

Tertiary Education Percentage Controls: ASTRA - 2021

Filtering CO2 outliers for calculating CO2 mean regressand

#Note: by removing CO2 emissions anomalies, I would need to assume that the anomalies are randomly distributed among 
#the dataset as I am still counting these vehicles in the number of registrations. This is unfortunately incorrect,
#as hybrids and certain brands appeared to be statistically more likely to be missing datapoints. Regarless, very 
#few vechiles are removed relative to the larger sample size so this impact is likely minimal

filteredNEUZU['RealTreibstoff_Code'] = filteredNEUZU['Treibstoff_Code']

print("Original dataset length: ",  "{:,}".format(len(filteredNEUZU)),"\n")

print("Filtering out values without energy category:", len(filteredNEUZU) - len(filteredNEUZU[filteredNEUZU['RealTreibstoff_Code'].notna()]), "terms")
CO2Filtered = filteredNEUZU[filteredNEUZU['RealTreibstoff_Code'].notna()]

print("Removing non-electric vehicles listed as having incorrectly low carbon emissions:",  len(CO2Filtered[((~CO2Filtered['RealTreibstoff_Code'].isin(['E','R','X'])) & (CO2Filtered['CO2'] <= 25))]), "terms")
CO2Filtered = (CO2Filtered[~((~CO2Filtered['RealTreibstoff_Code'].isin(['E','R','X'])) & (CO2Filtered['CO2'] <= 25))])

print("Replacing REX models with an engine as non-EVs:", len(CO2Filtered[(CO2Filtered['RealTreibstoff_Code'] == 'R') & (CO2Filtered['Hubraum'] > 0)]), "terms")
CO2Filtered['RealTreibstoff_Code'] = np.where((CO2Filtered['RealTreibstoff_Code'] == 'R') & (CO2Filtered['Hubraum'] > 0), 'REX-WithEngine', CO2Filtered['RealTreibstoff_Code'])

print("Setting CO2 emissions for EVs to zero:", len( (CO2Filtered[((CO2Filtered['RealTreibstoff_Code'].isin(['E','R','X'])) & (CO2Filtered['CO2'] > 10)) ] )), "terms")
CO2Filtered['CO2'] = np.where(CO2Filtered['RealTreibstoff_Code'].isin(['E','R','X']), 0, CO2Filtered['CO2'])

print("Filtering out values without CO2 emissions:",  "{:,}".format((len(CO2Filtered) - len(CO2Filtered[CO2Filtered['CO2'] >= 0]))), "terms") 
CO2Filtered = CO2Filtered[CO2Filtered['CO2'].notna()]

print("Filtering out incorrectly high CO2 emissions data:", len(CO2Filtered[(CO2Filtered['CO2'] > 500)]), "terms")
CO2Filtered = CO2Filtered[(CO2Filtered['CO2'] <= 500)]

print("\nCleaned CO2 data length:",  "{:,}".format(len(CO2Filtered)), " \n(which is a", round((len(filteredNEUZU) - len(CO2Filtered))/len(filteredNEUZU)*100,3), "percent decrease in size)" )

Original dataset length:  1,732,720 

Filtering out values without energy category: 0 terms
Removing non-electric vehicles listed as having incorrectly low carbon emissions: 199 terms
Replacing REX models with an engine as non-EVs: 156 terms
Setting CO2 emissions for EVs to zero: 770 terms
Filtering out values without CO2 emissions: 89,833 terms
Filtering out incorrectly high CO2 emissions data: 174 terms

Cleaned CO2 data length: 1,642,514  
(which is a 5.206 percent decrease in size)

Merging CO2 data with Regression Dataset

#Organizing CO2 emissions by canton year, quarter and merging with CO2 dataset
CO2Filtered['Year-Quarter'] = CO2Filtered['YearQuarter']

Co2Avg = CO2Filtered.groupby(['Inverkehrsetzung_Kanton','Year-Quarter']).CO2.agg(['mean', 'median', 'max', 'count']).reset_index().rename(columns={'mean':'CO2Mean', 'median':'CO2Median'})

Co2Avg['CantonYearQuarter'] = Co2Avg['Inverkehrsetzung_Kanton'].astype(str) + Co2Avg['Year-Quarter'].astype(str)
CYQ_preDummies_CO2 = pd.merge(left=CYQ_preDummies, right=Co2Avg, left_on='CantonYearQuarter', right_on='CantonYearQuarter').rename(columns={'CantonName':'Canton'})

Taking natural log of regressands

Note - this cell will run an error due to taking the natural log of 0s. The couple canton-year-quarter measurements without any EV registrations will be removed in the regressions

CYQ_preDummies_CO2['EVs_log']          = np.log(CYQ_preDummies_CO2['EVs'])
CYQ_preDummies_CO2['MkShr(100)_log']   = np.log(CYQ_preDummies_CO2['EVRatio'])
CYQ_preDummies_CO2['AllVehicles_log']  = np.log(CYQ_preDummies_CO2['AllVehicles'])
CYQ_preDummies_CO2['CO2Mean_log']      = np.log(CYQ_preDummies_CO2['CO2Mean'])

#removing natural logs of 0 in the regressions
CYQ_preDummies_CO2 = CYQ_preDummies_CO2.replace([np.inf, -np.inf], np.nan)
CYQ_preDummies_CO2 = CYQ_preDummies_CO2.drop(columns = ['Unnamed: 0', 'Year_y', 'Inverkehrsetzung_Kanton']).rename(columns = {"max": "Max CO2 emissions vehicle"})

C:\Users\lsvoi\anaconda3\lib\site-packages\pandas\core\arraylike.py:358: RuntimeWarning: divide by zero encountered in log
  result = getattr(ufunc, method)(*inputs, **kwargs)

# Vertical Data to csv file #
first_column = CYQ_preDummies_CO2.pop('CantonYearQuarter')
CYQ_preDummies_CO2.insert(0, 'CantonYearQuarter', first_column)
CYQ_preDummies_CO2.to_csv('EVData_VerticallyOrganized.csv')

Summary Statistic tables

regressands = ['EVs','EVRatio', 'AllVehicles', 'CO2Mean']
regressands_log = ['EVs_log','MkShr(100)_log', 'AllVehicles_log', 'CO2Mean_log']
controls = ['TaxRatio','Pct Tertiary', 'Population_Density_Over25', 'PopulationOver25']

#latex form commented out below
for tables in [regressands, regressands_log, controls]:
    display(  HTML(  informationalTable.correctLabels(CYQ_preDummies_CO2[tables].describe().round(2).to_html())  ))
    #print(CYQ_preDummies_CO2[tables].describe().round(2).to_latex().replace("mean","Mean").replace("count","Count").replace("std","St. Dev.").replace("min","Min").replace("max","Max").replace("624.00", "624").replace("622.00","622"))  

	Electric Vehicles	Pct. of Electric Vehicles	Total Number of Vehicles	Average CO2 Emissions
Count	624.00	624.00	624.00	624.00
Mean	124.82	4.77	2776.79	130.82
St. Dev.	222.57	4.93	2753.87	14.17
Min	0.00	0.00	176.00	73.17
25%	19.00	1.22	655.25	128.94
50%	47.00	2.60	1939.50	134.18
75%	133.00	6.68	3649.50	139.07
Max	2305.00	24.07	14052.00	153.22

	Electric Vehicles Registered, Logarithmic (Model 1)	Market Share of Electric Vehicles, Logarithmic (Model 2)	Total Vehicle Registrations, Logarithmic (Model 3)	Average Carbon Dioxide Emissions, Logarithmic (Model 4)
Count	622.00	622.00	624.00	624.00
Mean	3.84	1.02	7.42	4.87
St. Dev.	1.49	1.12	1.08	0.12
Min	0.00	-3.58	5.17	4.29
25%	2.94	0.20	6.48	4.86
50%	3.85	0.97	7.57	4.90
75%	4.89	1.90	8.20	4.94
Max	7.74	3.18	9.55	5.03

	Annual Taxes for EV vs. NonEV, Ratio	Tertiary Education Pct.	Population Den. (persons over 25/km)	Population (persons over 25)
Count	624.00	624.00	624.00	624.00
Mean	-0.53	32.54	359.86	237078.78
St. Dev.	0.60	6.66	756.40	250704.93
Min	-1.00	19.24	20.61	11034.00
25%	-1.00	27.99	65.30	52260.50
50%	-0.68	31.42	159.07	171714.50
75%	-0.42	35.93	253.37	292485.75
Max	1.92	52.01	3967.51	1143880.00

Verticallly Organized Dataset Sample

CYQ_preDummies_CO2[['GermanCantonName'] + ['YearQuarter'] + regressands + ['PP_1000'] + ['CS_1000'] + controls[:3]][372:565:4].rename(columns = {"EVs": "EVs Reg. (For Model 1)", "AllVehicles" : "All Vehicle Reg. (For Model 3) ", "CO2Mean": "Avg. Emission Level of Reg. (For Model 4)", "GermanCantonName": "Canton", "EVRatio" : "EV Registration % (For Model 2)", "PP_1000": "PP/1000", "CS_1000": "CS/1000", "TaxRatio": "Annual Tax, EV vs NonEV", "Population_Density_Over25": "Pop. Den. (25+)"})

	Canton	YearQuarter	EVs Reg. (For Model 1)	EV Registration % (For Model 2)	All Vehicle Reg. (For Model 3)	Avg. Emission Level of Reg. (For Model 4)	PP/1000	CS/1000	Annual Tax, EV vs NonEV	Pct Tertiary	Pop. Den. (25+)
372	Schaffhausen	2019-1	34	5.9028	576	137.455342	0.0	0.0	-0.125000	30.675235	203.573826
376	Schaffhausen	2020-1	25	5.7604	434	137.439716	0.0	0.0	-0.125000	29.856080	204.909396
380	Schaffhausen	2021-1	54	12.1896	443	117.584184	2.0	0.0	-0.125000	37.636936	207.043624
384	Schwyz	2016-1	24	1.4760	1626	140.332073	0.0	0.0	0.604706	30.063647	123.038546
388	Schwyz	2017-1	21	1.2419	1691	142.033313	0.0	0.0	0.604706	32.699502	124.936123
392	Schwyz	2018-1	39	2.5760	1514	145.528658	0.0	0.0	0.604706	32.769788	126.607930
396	Schwyz	2019-1	66	4.3854	1505	149.204068	0.0	0.0	0.604706	32.212421	128.530837
400	Schwyz	2020-1	77	6.0109	1281	144.883004	0.0	0.0	0.604706	33.036944	130.176211
404	Schwyz	2021-1	114	8.9062	1280	133.935943	0.0	0.0	0.604706	41.095971	132.296256
408	Solothurn	2016-1	25	1.1008	2271	131.318489	0.0	0.0	-1.000000	25.980065	246.984829
412	Solothurn	2017-1	24	1.0508	2284	133.529150	0.0	0.0	-1.000000	28.237685	250.236410
416	Solothurn	2018-1	63	2.8264	2229	132.872273	0.0	0.0	-1.000000	28.133925	252.786346
420	Solothurn	2019-1	98	4.5224	2167	137.574178	0.0	0.0	-1.000000	28.028625	255.116308
424	Solothurn	2020-1	83	4.7729	1739	134.964892	0.0	0.0	-1.000000	28.344835	257.586599
428	Solothurn	2021-1	132	7.8900	1673	121.940217	0.0	0.0	-1.000000	36.104096	260.367889
432	St. Gallen	2016-1	52	1.3259	3922	133.446170	0.0	0.0	-1.000000	26.348183	174.714215
436	St. Gallen	2017-1	65	1.5819	4109	134.003955	0.0	0.0	-1.000000	27.995172	176.995558
440	St. Gallen	2018-1	67	1.7521	3824	136.207958	0.0	0.0	-1.000000	27.955332	178.626357
444	St. Gallen	2019-1	175	4.6629	3753	140.209403	0.0	0.0	-1.000000	28.176141	180.167325
448	St. Gallen	2020-1	182	6.3393	2871	132.108194	0.0	0.0	-1.000000	29.404319	181.585884
452	St. Gallen	2021-1	310	10.2683	3019	118.023890	0.0	0.0	-1.000000	35.467715	183.335143
456	Thurgau	2016-1	33	1.4342	2301	132.548649	0.0	0.0	-0.882353	27.048487	193.511604
460	Thurgau	2017-1	41	1.7198	2384	134.084841	0.0	0.0	-0.882353	28.988833	197.175580
464	Thurgau	2018-1	35	1.5787	2217	135.164434	0.0	0.0	-0.882353	28.826809	200.225025
468	Thurgau	2019-1	166	7.1644	2317	134.604304	4.0	0.0	-0.882353	27.967584	202.915237
472	Thurgau	2020-1	148	8.6600	1709	130.805755	3.5	0.0	-0.882353	29.073131	205.719475
476	Thurgau	2021-1	213	12.0135	1773	115.251421	2.0	0.0	-0.882353	35.869201	208.733602
480	Ticino	2016-1	34	0.7532	4514	127.826155	0.0	0.0	-0.680000	30.154405	93.558321
484	Ticino	2017-1	31	0.6481	4783	130.173694	0.0	0.0	-0.680000	31.027426	94.366999
488	Ticino	2018-1	43	1.0007	4297	130.574399	0.0	0.0	-0.674015	32.443680	94.331081
492	Ticino	2019-1	90	2.0469	4397	137.297790	0.0	0.0	-0.674015	32.285259	94.437055
496	Ticino	2020-1	139	4.3114	3224	134.519129	2.0	0.5	-0.674015	31.973194	94.200213
500	Ticino	2021-1	179	4.9149	3642	124.998790	2.0	0.5	-0.674015	33.144302	94.248933
504	Uri	2016-1	2	0.6623	302	135.946488	0.0	0.0	-0.420973	19.236090	23.877437
508	Uri	2017-1	4	1.2987	308	135.117450	0.0	0.0	-0.420973	25.336991	23.996286
512	Uri	2018-1	1	0.3745	267	136.287356	0.0	0.0	-0.420973	22.642318	24.197772
516	Uri	2019-1	7	2.3729	295	143.771127	0.0	0.0	-0.420973	24.937951	24.331476
520	Uri	2020-1	7	3.2558	215	137.201878	0.0	0.0	-0.420973	22.245056	24.561746
524	Uri	2021-1	10	4.2735	234	126.095477	0.0	0.0	-0.420973	31.531350	24.729805
528	Valais	2016-1	25	0.7825	3195	139.880127	0.0	0.0	-0.600000	25.731333	46.315658
532	Valais	2017-1	51	1.6860	3025	139.109768	0.0	0.0	-0.600000	25.823717	46.965161
536	Valais	2018-1	24	0.8400	2857	141.920996	0.0	0.0	-0.600000	27.097610	47.531968
540	Valais	2019-1	85	3.0131	2821	147.247401	0.0	0.0	-0.600000	27.086451	48.063936
544	Valais	2020-1	86	3.8462	2236	142.287466	0.0	0.0	-0.600000	27.615261	48.407925
548	Valais	2021-1	212	8.9831	2360	124.903319	3.5	1.5	-0.600000	29.134491	49.017228
552	Vaud	2016-1	66	0.8852	7456	131.544263	0.0	0.0	-0.979625	37.310171	168.144147
556	Vaud	2017-1	94	1.2230	7686	132.074915	0.0	0.0	-0.979625	38.009681	171.161582
560	Vaud	2018-1	78	1.0338	7545	134.179463	0.0	0.0	-0.979625	38.990020	173.331880
564	Vaud	2019-1	294	3.8705	7596	136.786198	0.0	0.0	-0.979625	38.809152	174.961395

Converting dummy variables (at Canton and Year-Quarter level)

CYQ_dummies = pd.get_dummies(CYQ_preDummies_CO2, columns = ['YearQuarter'])
CYQ_dummies['Year-Quarter'] = CYQ_preDummies_CO2['YearQuarter']
CYQ_dummies = pd.get_dummies(CYQ_dummies, columns = ['Canton'])

Regression set-up

#lists for controls
priceControls =  ['CS_1000','TaxRatio']
socialControls = ['Pct Tertiary', 'PopDen(Over25)_log']
cantonControls = ["Canton_AG","Canton_AI","Canton_AR","Canton_BE","Canton_BL",
                  "Canton_BS","Canton_FR","Canton_GE","Canton_GL","Canton_GR",
                  "Canton_JU","Canton_LU","Canton_NE","Canton_NW","Canton_OW",
                  "Canton_SG","Canton_SH","Canton_SO","Canton_SZ","Canton_TG",
                  "Canton_TI","Canton_UR","Canton_VD","Canton_VS","Canton_ZG"]  # Canton_ZH excluded - dummy variable trap!
yearControls =   ["YearQuarter_2016-2","YearQuarter_2016-3","YearQuarter_2016-4",
                  "YearQuarter_2017-1","YearQuarter_2017-2","YearQuarter_2017-3","YearQuarter_2017-4",
                  "YearQuarter_2018-1","YearQuarter_2018-2","YearQuarter_2018-3","YearQuarter_2018-4",
                  "YearQuarter_2019-1","YearQuarter_2019-2","YearQuarter_2019-3","YearQuarter_2019-4",
                  "YearQuarter_2020-1","YearQuarter_2020-2","YearQuarter_2020-3","YearQuarter_2020-4",
                  "YearQuarter_2021-1","YearQuarter_2021-2","YearQuarter_2021-3","YearQuarter_2021-4"] #YearQuarter_2016-1 excluded - dummy variable trap!

#building the lists of the models
model1Regressors = ['PP_1000']
model2Regressors = ['PP_1000'] + priceControls
model3Regressors = ['PP_1000'] + priceControls + socialControls
model4Regressors = ['PP_1000'] + priceControls + socialControls + cantonControls
model5Regressors = ['PP_1000'] + priceControls + socialControls + yearControls
model6Regressors = ['PP_1000'] + priceControls + socialControls + cantonControls + yearControls

# Returns the OLS model with clustered StErrors based on given Xs/Y (regressors/regressand) inputs
def createModel(Xs,Y, cluster_id = -1):
    if type(cluster_id) == int:
        return sm.OLS(CYQ_dummies[[Y]], sm.add_constant(CYQ_dummies[Xs]), missing = 'drop').fit()    
    else: # clustered at canton level
        return sm.OLS(CYQ_dummies[[Y]], sm.add_constant(CYQ_dummies[Xs]), missing = 'drop').fit(cov_type='cluster', cov_kwds={'groups': cluster_id})
        
        
#Takes a regressand input and returns a summary table with regression output
def runRegression(regressand,i = 0):
    Y = regressand
    cluster_id =  CYQ_dummies.dropna(subset = [regressand])['GermanCantonName']
    
    #adding a constant term, dropping NaN ln(0) explained variables, clustering at canton level due to expected 
    #herteroskidasticity resulting from increase in variance with time) and using robust standard errros
    model1 = createModel(model1Regressors, Y, cluster_id)
    model2 = createModel(model2Regressors, Y, cluster_id)
    model3 = createModel(model3Regressors, Y, cluster_id)
    model4 = createModel(model4Regressors, Y, cluster_id)
    model5 = createModel(model5Regressors, Y, cluster_id)
    model6 = createModel(model6Regressors, Y, cluster_id)
    
    dfmodels = summary_col([model1, model2,model3,model4,model5,model6],
        model_names=[f'\nPurchase Premium     \n({i}.1)         \n',
                     f'\nVehicle Price Controls    \n({i}.2)    \n',
                     f'\nSocial Controls     \n({i}.3)          \n',
                     f'\nCanton FE        \n({i}.4)             \n',
                     f'\nTime FEs w/out Canton FEs    \n({i}.5) \n',
                     f'\nTime FEs & Canton FEs\n({i}.6)         \n'], 
        info_dict={
            'Number of Observations':lambda x: "{0:d}".format(int(x.nobs)),
            'Canton Fe': lambda x: {1: 'No', 3: 'No', 5: 'No', 30: 'Yes', 28: 'No', 53: 'Yes'}[int(x.df_model)],
            'Time Fe':  lambda x: {1: 'No', 3: 'No', 5: 'No', 30: 'No', 28: 'Yes', 53: 'Yes'}[int(x.df_model)]},
        float_format='%0.3f', stars=True,
        regressor_order=model3Regressors, 
        drop_omitted=True)   

    return dfmodels

Running Regressions

1) The natural log of electric vehicle sales by quarter-year-canton

2) The natural log of electric vehicle market share by quarter-year-canton

3) The natural log of all vehicle sales by quarter-year-canton

4) The natural log of the average CO2 emissions of newly registered vehicles within each quarter-year-canton

TypeLog_Label = {'EVs_log': 'Electric Vehicles Registered, Logarithmic','MkShr(100)_log':'Market Share of Electric Vehicles, Logarithmic','AllVehicles_log':'Total Vehicle Registrations, Logarithmic', 'CO2Mean_log':'Average Carbon Dioxide Emissions, Logarithmic'}
modelNum = 0
for regressand in regressands_log:
    modelNum += 1
    #print(regressand, runRegression(regressand).as_latex().replace('***   ','^{***}').replace('**   ','^{**}').replace('*   ','^{*}').replace("R-squared",'R\ '+ 'textsuperscript{2}').replace("N",'R\ '+ 'textsuperscript{2}').replace('& 4,6               & 4,6                    & 4,6                       & 4,6                        & 4,6                       & 4,6                     \\', '& No               & No                    & No                       & Yes                       & No                        & Yes                     \\').replace('& 5,6               & 5,6                    & 5,6                       & 5,6                        & 5,6                       & 5,6                     \\', '& No               & No                    & No                       & No                        & Yes                       & Yes                     \\').replace('\_', ' /'))
    print(f' Model {modelNum}:\n', TypeLog_Label[regressand])
    display(HTML(runRegression(regressand,modelNum).as_html().replace("_", " / ").replace("PopDen(Over25) / log","Ln of Pop. Den.").replace("Pct Tertiary", "Tertiary Education Pct.")))
    print("\n", "\n", "\n")

 Model 1:
 Electric Vehicles Registered, Logarithmic

	Purchase Premium	Vehicle Price Controls	Social Controls	Canton FE	Time FEs w/out Canton FEs	Time FEs & Canton FEs
	(1.1)	(1.2)	(1.3)	(1.4)	(1.5)	(1.6)
PP / 1000	0.527***	0.346***	0.493***	0.186**	0.285***	0.086***
	(0.111)	(0.100)	(0.082)	(0.092)	(0.094)	(0.016)
CS / 1000		1.214***	0.559**	0.335***	0.308	0.038
		(0.231)	(0.275)	(0.124)	(0.330)	(0.051)
TaxRatio		-0.274	-0.307	-0.281**	-0.305	-0.103
		(0.323)	(0.290)	(0.133)	(0.271)	(0.063)
Tertiary Education Pct.			0.119***	0.093***	0.064	0.005
			(0.034)	(0.026)	(0.043)	(0.018)
Ln of Pop. Den.			-0.018	27.502***	0.217	-4.800
			(0.243)	(4.144)	(0.297)	(4.008)
R-squared	0.047	0.101	0.361	0.862	0.514	0.950
R-squared Adj.	0.046	0.097	0.356	0.855	0.491	0.945
Number of Observations	622	622	622	622	622	622
Canton Fe	No	No	No	Yes	No	Yes
Time Fe	No	No	No	No	Yes	Yes

 
 

 Model 2:
 Market Share of Electric Vehicles, Logarithmic

	Purchase Premium	Vehicle Price Controls	Social Controls	Canton FE	Time FEs w/out Canton FEs	Time FEs & Canton FEs
	(2.1)	(2.2)	(2.3)	(2.4)	(2.5)	(2.6)
PP / 1000	0.450***	0.360***	0.455***	0.202*	0.164***	0.099***
	(0.053)	(0.057)	(0.043)	(0.106)	(0.056)	(0.029)
CS / 1000		0.612***	0.169	0.386**	-0.195*	0.107
		(0.145)	(0.172)	(0.158)	(0.104)	(0.106)
TaxRatio		-0.121	-0.068	-0.298**	-0.073	-0.101*
		(0.134)	(0.156)	(0.127)	(0.120)	(0.055)
Tertiary Education Pct.			0.106***	0.113***	0.022	0.029
			(0.021)	(0.035)	(0.016)	(0.034)
Ln of Pop. Den.			-0.287***	29.928***	0.077	-5.518
			(0.091)	(6.818)	(0.070)	(4.371)
R-squared	0.061	0.084	0.284	0.701	0.749	0.865
R-squared Adj.	0.060	0.080	0.278	0.686	0.737	0.852
Number of Observations	622	622	622	622	622	622
Canton Fe	No	No	No	Yes	No	Yes
Time Fe	No	No	No	No	Yes	Yes

 
 

 Model 3:
 Total Vehicle Registrations, Logarithmic

	Purchase Premium	Vehicle Price Controls	Social Controls	Canton FE	Time FEs w/out Canton FEs	Time FEs & Canton FEs
	(3.1)	(3.2)	(3.3)	(3.4)	(3.5)	(3.6)
PP / 1000	0.078	-0.013	0.040	-0.016	0.121	-0.013
	(0.094)	(0.086)	(0.091)	(0.020)	(0.102)	(0.019)
CS / 1000		0.606***	0.390	-0.049	0.503*	-0.068
		(0.174)	(0.264)	(0.051)	(0.262)	(0.070)
TaxRatio		-0.150	-0.236	0.014	-0.228	-0.003
		(0.268)	(0.229)	(0.019)	(0.240)	(0.043)
Tertiary Education Pct.			0.013	-0.020*	0.043	-0.024
			(0.026)	(0.011)	(0.040)	(0.019)
Ln of Pop. Den.			0.268	-2.496	0.140	0.766
			(0.230)	(3.194)	(0.263)	(1.267)
R-squared	0.002	0.029	0.143	0.947	0.193	0.960
R-squared Adj.	0.000	0.024	0.136	0.945	0.155	0.956
Number of Observations	624	624	624	624	624	624
Canton Fe	No	No	No	Yes	No	Yes
Time Fe	No	No	No	No	Yes	Yes

 
 

 Model 4:
 Average Carbon Dioxide Emissions, Logarithmic

	Purchase Premium	Vehicle Price Controls	Social Controls	Canton FE	Time FEs w/out Canton FEs	Time FEs & Canton FEs
	(4.1)	(4.2)	(4.3)	(4.4)	(4.5)	(4.6)
PP / 1000	-0.042***	-0.030***	-0.035***	-0.017***	-0.008**	-0.008**
	(0.013)	(0.011)	(0.009)	(0.006)	(0.004)	(0.003)
CS / 1000		-0.084***	-0.057*	-0.100***	-0.003	-0.002
		(0.026)	(0.032)	(0.011)	(0.013)	(0.005)
TaxRatio		0.009	0.004	-0.009	0.016**	0.013***
		(0.009)	(0.013)	(0.006)	(0.007)	(0.004)
Tertiary Education Pct.			-0.007**	-0.027***	0.005***	0.000
			(0.003)	(0.003)	(0.002)	(0.001)
Ln of Pop. Den.			0.025*	0.076	-0.029***	-1.261***
			(0.013)	(0.409)	(0.008)	(0.443)
R-squared	0.044	0.074	0.141	0.519	0.905	0.951
R-squared Adj.	0.042	0.070	0.134	0.494	0.901	0.946
Number of Observations	624	624	624	624	624	624
Canton Fe	No	No	No	Yes	No	Yes
Time Fe	No	No	No	No	Yes	Yes

Making some graphs

# Makes new folders that store the following images as png files
countryLevelGraphs.createNewFolders(cantons_abbr)

Average CO2 emissions vs. EV Market Share

countryLevelGraphs.CO2vsEVMarketShare(CO2Filtered)

Total Number of EVs, Non EVs, EVs

countryLevelGraphs.EVsNonEVsAllVehicles()

Plotting Electric Vehicle Market Share by Canton over Time

for c in ['SZ','ZH','TG','TI', 'SH', 'VS']:
    cantonLevelGraphs.plotEVRatio(c)

Plotting by number of EV sales

for c in ['SZ','ZH','TG','TI', 'SH', 'VS']:
    cantonLevelGraphs.plotEVs(c)

Plotting by total number of vehicle sales

for c in ['SZ','ZH']:
    cantonLevelGraphs.allVehicles(c)

Plotting the residuals at the canton level for a particular control-group model

This residual analysis allows for better visualization of regression output meanings

### creates graph for analysis based on a given canton abbreviation input ###
def stage12346(c_abbr, Y = 'MkShr(100)_log'): 
    residGraph.graph("","",                                                                                                      c_abbr, Y, 0, True )
    residGraph.graph( createModel(model1Regressors[1:], Y),  "(controlled by a constant)",                                       c_abbr, Y, 1, False)
    residGraph.graph( createModel(model2Regressors[1:], Y),  "(controlled by a constant, vehicle price)",                        c_abbr, Y, 2, False)
    residGraph.graph( createModel(model3Regressors[1:], Y),  "(controlled by a constant, vehicle price, social)",                c_abbr, Y, 3, False)
    residGraph.graph( createModel(model4Regressors[1:], Y),  "(controlled by a constant, vehicle price, social, canton)",        c_abbr, Y, 4, False)
    residGraph.graph( createModel(model6Regressors[1:], Y),  "(controlled by a constant, vehicle price, social, canton, quarter)",  c_abbr, Y, 6, False)
    
    ## Thurgau graph well explains variation post-policy implimentation, however Schauffhausen explains this variation less so
    ## These can be seen when comparing the full model .61 without purchase premium and .62 with purchase premium control.
    if c_abbr in ['TG', 'SH']:
        residGraph.graph(createModel(model6Regressors[1:], Y),"(controlled by a constant, vehicle price, social, canton, quarter)", c_abbr, Y, 61, True)
        residGraph.graph(createModel(model6Regressors, Y),    "(controlled by all controls and the purchase premiums policies)", c_abbr, Y, 62, True)

### builds photos into a gif ###
def makeGifs(cantonName):
    images = []
    for i in [1,2,3,4,6,6,6]:
        images.append(imageio.imread("2) Gifs/" + cantonName + '/img' + cantonName + str(i) + '.png'))
    exportname = cantonName +'.gif'
    imageio.mimsave("2) Gifs/" + cantonName + '/' + exportname, images, duration = 1.9)
    
### stages and creates the gif ###
def stageAndGifs(cantonName):
    stage12346(cantonName)
    makeGifs(cantonName)

for cantonName in ['TG', 'SH', 'VS', 'TI']:
    stageAndGifs(cantonName)

# Display some of the gifs
def displayGif(abbr):
    display(Image(data=open('2) Gifs/' + abbr + '/' + abbr + '.gif','rb').read(), format='png'))
for abbr in [ 'TG', 'TI', 'SH', 'VS']:
    displayGif(abbr)

print("Unexplained residuals indicate positive effect of purchase premium policy in canton of Valais, \nshown by green region within the red circle in image below")
Image(filename="2) Gifs/VS/imgVS6.png")

Unexplained residuals indicate positive effect of purchase premium policy in canton of Valais, 
shown by green region within the red circle in image below

### builds photos into a gif ###
def prePostResiduals(cantonNames):
    for cantonName in cantonNames:
        images = []
        for i in [61,62]:
            images.append(imageio.imread("2) Gifs/" + cantonName + '/img' + cantonName + str(i) + '.png'))
        imageio.mimsave("2) Gifs/" + cantonName + '/' + cantonName +'ResidualChange_PrePost.gif', images, duration = 1.2)

prePostResiduals(['TG'])
print("Flattening of residuals post-policy implimentation indicates well-explained residuals and a positive effect of purchase premium policy in canton of Thurgau")

for abbr in ['TG']:
    display(Image(data=open('2) Gifs/' + abbr + '/' + abbr + 'ResidualChange_PrePost.gif','rb').read(), format='png'))

Flattening of residuals post-policy implimentation indicates well-explained residuals and a positive effect of purchase premium policy in canton of Thurgau

Thank you for your time.

-Lucas