eststo clear


if("`c(os)'"=="Windows"){
	local tableSaveDir = "Tables\DiD\"
	local plotSaveDir = "Plots\"
	local pythonSaveDir = "PythonScripts\"
}
else{
	local tableSaveDir = "Tables/DiD/"
	local plotSaveDir = "Plots/"
	local pythonSaveDir = "PythonScripts/"
}

//Establish Control Variable Sets - VMT
local Spec1ControlsVMT = "ib51.stateGroup ib1995.year c.normYear#i.stateGroup asinhMetroPop asinhNonMetroPop asinhMetroInc asinhNonMetroInc logrealmeangasprice logemployment loglicenseddrivers logrealstategdp logroadmileage"
local Spec2ControlsVMT = "`Spec1ControlsVMT'"
local Spec3ControlsVMT = "`Spec2ControlsVMT' i.vmtDummyExtSpec3"
local Spec4ControlsVMT = "`Spec2ControlsVMT' i.vmtDummyExtSpec4"


//Establish Control Variable Sets - Gas
local Spec1ControlsGas = "ib51.stateGroup ib1995.year c.normYear#i.stateGroup asinhMetroPop asinhNonMetroPop asinhMetroInc asinhNonMetroInc logrealmeangasprice logemployment loglicenseddrivers logrealstategdp logroadmileage"
local Spec2ControlsGas = "`Spec1ControlsGas'"
local Spec3ControlsGas = "`Spec2ControlsGas' i.gasUseDummyExtSpec3"
local Spec4ControlsGas = "`Spec2ControlsGas' i.gasUseDummyExtSpec4"


//Cluster variable declaration
local clusterVar = "stateGroup"

//Treatment variable 
local treatmentVar = "1.nosafetyind"

//Run Spec2 regressions - prep for suest command

qui regress logvmt i.nosafetyind `Spec2ControlsVMT'
scalar travelK = e(df_m)+1
estimates store travelEq


qui regress loghighwaygasuse i.nosafetyind `Spec2ControlsGas'
scalar gasK = e(df_m)+1
estimates store gasEq

//Simultaneously estimate reduced form esquations
qui suest travelEq gasEq, vce(cluster `clusterVar')
scalar nClust = e(N_clust)
scalar N = e(N)
matlist e(b)["y1", "travelEq_mean:`treatmentVar'"]
scalar delta = e(b)["y1", "travelEq_mean:`treatmentVar'"]+0
scalar theta = e(b)["y1", "gasEq_mean:`treatmentVar'"]
scalar deltaVarUncorrected = e(V)["travelEq_mean:`treatmentVar'", "travelEq_mean:`treatmentVar'"]
scalar thetaVarUncorrected = e(V)["gasEq_mean:`treatmentVar'", "gasEq_mean:`treatmentVar'"]
scalar covarUncorrected = e(V)["gasEq_mean:`treatmentVar'", "travelEq_mean:`treatmentVar'"]

//Collect all necessary statistics
scalar travelCorrection = (N-1)/(N-travelK+(nClust-1))
scalar gasCorrection = (N-1)/(N-gasK+(nClust-1))
scalar deltaVar = deltaVarUncorrected*travelCorrection
scalar thetaVar = thetaVarUncorrected*gasCorrection
scalar covarHigh = max(travelCorrection*covarUncorrected,gasCorrection*covarUncorrected)
scalar covarLow = min(travelCorrection*covarUncorrected,gasCorrection*covarUncorrected)
scalar covar = covarLow

//Calculate wald-test stastics
matrix r = [1, -1]

matrix betaMat = [delta \ theta]
matrix covMat = [deltaVar, covar \ covar, thetaVar]

matrix outer = r * betaMat
matrix inner = r * covMat * r'

matrix chi2Val = outer*inv(inner)*outer'

scalar pVal = 1-chi2(1, chi2Val[1, 1])
//Display p-value
di pVal
