/*  1. Forecast FFT model

     In this do file we prepare the dataset used to estimate our scheme level FFT models.
	 This do file runs the regression for FFT enhancement totex model (forecast data model using totex and capacity increase). 
     Each FFT 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\FFT\Forecast"

import excel "Wastewater - Storm overflows; enhancement expenditure model.xlsx", sheet("Stata dataset (real) - FFT For") cellrange(A1:H104) 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 source 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 source                             /// Panel data identifiers 
			l_s_increase   	                                        /// Size of capacity increase

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

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

    * Import data
	
    use "FFT 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

	* Drop instances of missing or zero capacity increase/totex (if any)
	
	drop if l_s_increase == 0
	drop if missing(l_s_increase)
	
	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 source temp1
	
	* Transform to natural logs
	
	foreach x of varlist temp1-temp2 {
		g ln`x' = ln(`x')
	} 

	drop temp* lntemp*  

 *=================================================================
 *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 l_s_increase 														"liters per second increase"
	

 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			l_s_increase				    
   global reg2   lntotex        lnl_s_increase
   
	* 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}
}
	
* Predict Cook's distances and drop all observations with a Cook's distance >= 4 / N 
 
     reg totex l_s_increase
	   predict cook1, cooksd, if e(sample)  
	  list companycode stw totex if cook1 >= 4 / _N  
	  
	 reg lntotex lnl_s_increase
	   predict cook2, cooksd, if e(sample)  
	  list companycode stw totex if cook2 >= 4 / _N  
	  
	  g outlier = 0
	  replace outlier = 1 if cook1 >= 4 / _N  | cook2 >= 4 / _N 
	  
* Export to excel the full dataset to also keep cost and cost driver information of cooks distance outliers
	
	export excel companycode stw soid  totex l_s_increase outlier using "FFT - Costs and costs drivers.xlsx", sheet("Costs and costs drivers") sheetmodify firstrow(variables)
	
 *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
	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 "FFT - 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 "FFT - coefficients datestamp.xls", cell (A1) sheetmodify firstrow(variables)
	
	*Export the results
		
	estout using "FFT - 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)				 
		
		
*************************************************************************************************************