/*  1. Forecast STW storm overflow model

     In this do file we prepare the dataset used to estimate our scheme level SO STW models.
	 This do file runs the regression for SO enhancement totex model (forecast data model using totex and total storage). 
     Each SO totex model is run in a different do file to ensure the correct Cook's distance outliers of the relevant model are excluded and to reflect that there are two datasets (forecast and historical).
   
*/

*This do file has the following stages:

*1.- Data preparation - Import the master dataset from Excel and assign labels and codes for all variables. 
 
*2.- Variables generation using codes 

*3. Structure into a panel dataset format

*4.- Data cleaning

*5.- prepare macros for regressions

*6.- Run regressions and tests

   *6a) OLS - models for outlier exclusion
   *6b) OLS - models

*7.- Export results and matrices into Excel


 *======================================================================================== 
 *1.- Data preparation
 *========================================================================================

   *1a.- Import the master dataset from Excel and assign labels and codes for all variables. 
clear
matrix drop _all
set more off
set matsize 800
estimates drop _all
macro drop _all

cd "O:\OFWSHARE\Cost assessment\PR24\PR24 Final Determinations\Wastewater enhancement\Final Determinations Run\Storm overflows\STW\Forecast"

import excel "Wastewater - Storm overflows; enhancement expenditure model.xlsx", sheet("Stata dataset (real) - SO For") cellrange(A1:M2530) firstrow

	*Next, we assign names and labels with the first and second rows respectively to each variable using the following loop: 
	
foreach var of varlist * {
  label  variable  `var' "`=`var'[1]'"
}

drop in 1
	
	*Then can destring all variables 
	 *We generate two variables, "temp1" and "temp2" so we always destring the relevant variables no matter what order we have:
g temp1		= ""
g temp2		= ""
order 		companycode soid stw site_type model_type completion_date temp1
destring 	temp1-temp2, replace force
drop 		temp1 temp2
drop if soid=="0"


    
	

*We check the number of osbervations per company and financial year
		
tab companycode, m	

 *========================================================================================
 *2.- Variables generation using the codes  
 *======================================================================================== 
   

  	* The next step is to keep the relevant variables for logs and real terms transformation where applicable.	
	
	keep	opex capex												/// Cost variables in real terms
		    companycode soid stw site_type model_type completion_date              /// Panel data identifiers 
			total_storage grey_storage green_storage other_storage 	              /// Size of storage

			
  ********************************************************************************************************************************************************************************
  *2.c - SAVE DATASET
	save "SO APP20 dataset", replace 

 
 *===================================================
 *3. Structure into a panel dataset format
 *===================================================

    * Import data
	
    use "SO APP20 dataset", clear

	* Assign numerical values to company codes and eaids:
		
    encode companycode, g(companycode_)
    encode soid, g(soid_)
	
 *===================================================
 *4. Data cleaning
 *===================================================

	* We will run the model on SO enhancement totex.
	
	gen totex = capex + opex
	
	* For now we only want to run network models so we drop STWs
	
	drop if site_type == "Network"
	
	* We want to run our models on Network Grey/Hybrid models and in this instance for the purpose of comparison we will also include the STWs which have a FFT and storage driver. 
	
	keep if model_type	== "Grey/Hybrid" | model_type == "FFT/GreyHybrid"
	
	* Replace total storage units
	
	replace total_storage = total_storage / 100
	
	* Drop instances of missing or zero storage/totex (if any)
	
	drop if total_storage == 0
	drop if missing(total_storage)
	
	drop if totex <= 0
	drop if missing(totex)
	
    * Transform variables to log. We don't transform those variables expressed in percentage. 
    
	g temp1=1
	g temp2=1
	
	* Then place those variables in the left that are not being transformed 
	
	order companycode soid stw site_type model_type completion_date temp1
	
	* Transform to natural logs
	
	foreach x of varlist temp1-temp2 {
		g ln`x' = ln(`x')
	} 
	
	drop temp* lntemp*

	*Need to convert date variable to a stata date format
	gen completion_date_s = date(completion_date, "MDY")
	
	*Changing numerical date into standard calendar date format
	format completion_date_s %td
	
	* Extracting year from date
	gen completion_year = year(completion_date_s)
	
	* Creating date dummy variable
	gen year_dummy = (completion_year > 2030)

    * Generate an outlier dummy to populate with Cook's distance outliers
	gen outlier = 0

 *=================================================================
 *5.- Prepare the macros for different regression specifications 
 *=================================================================  

 * Label variables for easy understanding
	
	* Dependent variables
	label var opex 														        "Opex, £m (2022/23 prices)" 
	label var capex                                                             "Capex, £m (2022/23 prices)" 
	label var totex                                                             "Totex, £m (2022/23 prices)"
	
	* Panel variables
	label var companycode 														"Company code"
	label var stw																"Name of the sewage treatment work"
									
	* Independent variables
	label var total_storage 														"Total storage, meters cube"
	

 global regressionList = "1 2"	// set the list of regressions to run, if you do not want to run a particular regression then take the number out.

	* Models
	
   global reg1   totex			total_storage				    
   global reg2   lntotex        lntotal_storage
   
	* Create a global list of all companies to calculate efficiency scores
	global companyList =  "ANH NES NWT SRN SVH SWB TMS WSH WSX YKY"  				// do not change

 *=================================================================	
 *6.- Models' analysis
 *=================================================================
 *6a) OLS - models for outlier exclusion
foreach num of numlist $regressionList {
	global i=`num'

	*Run Pooled OLS with clustered standard errors at company level
	eststo  ols${i}: reg ${reg$i}
	estadd scalar RSS = e(rss): ols${i}
	local RSS = e(rss)
	local df = e(df_m)
	predict residuals${i},resid
		
	*Adjusted R squared
	estadd scalar R_squared_adj = e(r2_a): ols${i}
	estadd scalar R_squared = e(r2): ols${i}

	*VIF test
	vif
	estadd scalar VIF_statistic = r(vif_1): ols${i}	
 
	*RESET Test
	ovtest
	estadd scalar RESET_P_value = r(p): ols${i} 
	
	*Normality Test
	sktest residuals${i}
	estadd scalar Normality = r(P_chi2): ols${i}
	
		
	* Heteroscadsticity test
	quietly eststo olshet${i}: reg ${reg$i}
	estat hettest 
	estadd scalar Heteroscedasticity = r(p): ols${i}
	eststo drop olshet${i}	

	*Add title to the regression
	estadd local Econometric_model = "Pooled OLS": ols${i}
	
	* Drop calculated model variables, to save on memory
	drop _est_ols${i} residuals${i}
}

save "SO Data w AMP9 Schemes", replace
	
* Predict Cook's distances and drop all observations with a Cook's distance >= 4 / N 

     tab year_dummy

     drop if year_dummy == 1
 
     reg totex total_storage 
	   predict cook1, cooksd, if e(sample)  
	  list companycode stw totex if cook1 >= 4 / _N  
	  
	 reg lntotex lntotal_storage 
	   predict cook2, cooksd, if e(sample)  
	  list companycode stw totex if cook2 >= 4 / _N  
	  

	  replace outlier = 1 if cook1 >= 4 / _N  | cook2 >= 4 / _N 					

* Export to excel using the new set of outliers

	export excel companycode stw soid  totex total_storage outlier year_dummy completion_date_s using "SO - Costs and costs drivers.xlsx", sheet("Costs and costs drivers") sheetmodify firstrow(variables)
	
	clear
* Import the corect set of outliers, then merge with the full dataset (including the AMP9 sites)	
	import excel using "SO - Costs and costs drivers.xlsx", firstrow
	
	keep soid outlier
	
	merge 1:1 soid using "SO Data w AMP9 Schemes"
	
	tab year_dummy
	
* Export to excel the full dataset to also keep cost and cost driver information of cooks distance outliers
	
	export excel companycode stw soid  totex total_storage outlier year_dummy completion_date_s using "SO - Costs and costs drivers.xlsx", sheet("Costs and costs drivers") sheetmodify firstrow(variables)

 *6b) OLS - models
	
 *6b) OLS - models

foreach num of numlist $regressionList {
	global i=`num'

	*Run Pooled OLS with clustered standard errors at company level
	eststo  ols${i}: reg ${reg$i} if outlier == 0 & year_dummy == 0
	estadd scalar RSS = e(rss): ols${i}
	local RSS = e(rss)
	local df = e(df_m)
	predict residuals${i},resid

	*Adjusted R squared
	estadd scalar R_squared_adj = e(r2_a): ols${i}
	estadd scalar R_squared = e(r2): ols${i}

	*VIF test
	vif
	estadd scalar VIF_statistic = r(vif_1): ols${i}	
 
	*RESET Test
	ovtest
	estadd scalar RESET_P_value = r(p): ols${i} 
	
	*Normality Test
	sktest residuals${i}
	estadd scalar Normality = r(P_chi2): ols${i}
	
		
	* Heteroscadsticity test
	quietly eststo olshet${i}: reg ${reg$i}
	estat hettest 
	estadd scalar Heteroscedasticity = r(p): ols${i}
	eststo drop olshet${i}	

	*Add title to the regression
	estadd local Econometric_model = "Pooled OLS": ols${i}
	
	* Drop calculated model variables, to save on memory
	drop _est_ols${i} residuals${i}
}
	
*==================================================
 *7.- Export results and matrices into Excel	
 *==================================================
 	*Export a matrix for the coefficients that are used in feeder models 2 and 4.
	estout * 		using "SO - coefficients.xls"	, replace ///
					stats(depvar)


	g run_time_stamp = "$S_TIME" 
	g run_date_stamp = "$S_DATE"
	
	export excel run_date_stamp run_time_stamp using "SO - coefficients datestamp.xls", cell (A1) sheetmodify firstrow(variables)
	
	*Export the results
		
	estout using "SO - Results.xls", replace 														///
		cells(b(star fmt(3)) p(par({ }))) starlevels( * 0.10 ** 0.05 *** 0.010) 						///
		stats(Econometric_model depvar N vce R_squared R_squared_adj RESET_P_value VIF_statistic Normality Heteroscedasticity) mlabels(,titles)				 
		
		
*************************************************************************************************************