/* 3. Historical model - sq enhanced consent v2.0 In this do file we prepare the dataset used to estimate our scheme level phosphorus removal models. This do file runs the regression for phosphorus enhancement totex model PR3 (historical data model using PE, historical consent, enhanced consent and enhanced consent squared). Each phosphorus enhancement 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\Phosphorus\3. Historical model - squared enhanced consent" import excel "PR24CA60 - WW - P-removal.xlsx", sheet("Stata dataset (real) historical") cellrange(A1:L763) 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 eaid stw completion_date temp1 destring temp1-temp2, replace force drop temp1 temp2 drop if eaid=="0" *We check the number of osbervations per company and financial year tab companycode, m *======================================================================================== *2.- Variables generation using the codes *======================================================================================== * Denote the scheme completion year. gen completion_year = real(substr(completion_date, strlen("`string'")-4, 4)) * 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 eaid stw completion_year /// Panel data identifiers and completion year pe_served pe_design /// Size of works variables historical_consent enhanced_consent /// Treatment complexity variables transfer /// Dummy for transfers ******************************************************************************************************************************************************************************** *2.c - SAVE DATASET save "P 7F dataset", replace *=================================================== *3. Structure into a panel dataset format *=================================================== * Import data use "P 7F dataset", clear * Assign numerical values to company codes and eaids: encode companycode, g(companycode_) encode eaid, g(eaid_) *=================================================== *4. Data cleaning *=================================================== * We will run the model on P-removal enhancement totex. gen totex = capex + opex * We drop observations where enhanced consents are equal to zero drop if enhanced_consent == 0 * We assume that STWs that do not treat phosphorus can be considered to have consents of 5mg/l. g consent_assumed_historical = historical_consent replace consent_assumed_historical = 5 if historical_consent == 0 * We also generate a consent difference variable. Consent change variable is dropped where it is higher than zero (ie there is no tightened consent) g consent_change = consent_assumed_historical - enhanced_consent drop if consent_change <= 0 * Drop transferring from works drop if transfer == 1 * Drop PE that are equal to zero or missing drop if pe_served == 0 drop if pe_served == . * Drop if totex is zero drop if totex == 0 * Replace PE variables with thousands replace pe_served = pe_served / 1000 replace pe_design = pe_design / 1000 * Generate an outlier dummy to populate with Cook's distance outliers gen outlier = 0 * Dummy variable denoting if a STW has P consent <= 0.25 mg/l (TAL dummy) g p_below_025mgl = 0 replace p_below_025mgl = 1 if enhanced_consent <= 0.25 * Calculate square of P consent for non-linear modelling g sq_enhanced_consent = enhanced_consent^2 * Dummy variable denoting if a STW has P consent >= 2 mg/l (optimisations) g enhanced_consent_gr2 = 0 replace enhanced_consent_gr2 = 1 if enhanced_consent >= 2 tab enhanced_consent_gr2 * Drop if a STW has P consent >= 2 mg/l (optimisations) drop if enhanced_consent_gr2 == 1 * Export to excel export excel companycode stw eaid completion_year totex pe_served consent_assumed_historical enhanced_consent sq_enhanced_consent p_below_025mgl outlier enhanced_consent_gr2 using "P - Costs and costs drivers.xlsx", sheet("Costs and costs drivers") sheetmodify firstrow(variables) *================================================================= *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 pe_served "Population equivalent served, 000s" label var pe_design "Design population equivalent, 000s" label var historical_consent "Historical Phosphorus consent, mg/l" label var consent_assumed_historical "Historical Phosphorus consent with no consent replaced with 5mg/l, mg/l" label var enhanced_consent "Enhanced Phosphorus consent, mg/l" label var sq_enhanced_consent "Enhanced Phosphorus consent squared, mg/l" label var p_below_025mgl "Dummy variable denoting if Enhanced Phosphorus consent <= 0.25mg/l (TAL dummy)" global regressionList = "1" // 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 pe_served consent_assumed_historical enhanced_consent sq_enhanced_consent * 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 - predict Cook's distances and drop all observations with a Cook's distance >= 4 / N reg totex pe_served consent_assumed_historical enhanced_consent sq_enhanced_consent predict cook, cooksd, if e(sample) list companycode stw eaid totex if cook >= 4 / _N replace outlier = 1 if cook >= 4 / _N * Export to excel using the new set of outliers export excel companycode stw eaid completion_year totex pe_served consent_assumed_historical enhanced_consent sq_enhanced_consent p_below_025mgl outlier enhanced_consent_gr2 using "P - 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 excluding outliers 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 estout * using "P - 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 "P - coefficients datestamp.xls", cell (A1) sheetmodify firstrow(variables) *Export the results estout using "P - 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) *************************************************************************************************************