#rm(list = ls())
setwd(dirname(rstudioapi::getActiveDocumentContext()$path)) 
#set directory to source file
list.files()
library(spatialreg)
library(spdep)
#sourcing files for data and functions - takes a bit of time
source("Create_Dat3_final.R")
source("Funciones.R")
#----------------------------------------------------------------------------------
names(datk)
#-------------------------------------------
# Function to run sub-samples
# Arguments:  var.sub - the discrete variable along which sub-samples are created
#             formula - the formula to run on the sub-sample; defaults to fml.1
# This function only prints out output

ssregs=function(var.sub,formula=fml.1){
  u.var = unique(var.sub)
  for (l in 1:length(u.var)) {
    cat("Sub-sample is ",u.var[l]," \n")
    id = which(var.sub==u.var[l])
    Modl = regSARAR(formula = formula,subsamp = id) #use sub-sample
    summary.g2sls(Modl,marg_effect=FALSE,bet.id=3,L=1000,cl=5)
    cat(rep("=",45)); cat("\n")
  }#end for l
}

#=====================================================================================>
#********* Begin Printing out output to file
Out_log <- file("Out_log3_final_SARAR.txt") # File name of output log

sink(Out_log, append = TRUE, type = "output") # Writing console output to log file
sink(Out_log, append = TRUE, type = "message")

cat(readChar(rstudioapi::getSourceEditorContext()$path, # Writing currently opened R script to file
             file.info(rstudioapi::getSourceEditorContext()$path)$size))
#=====================================================================================>



#----------------------------------------------------------------------------------
cat("\n")
cat("Full Sample - Results","\n")
fml.1 = as.formula(paste("wrk~pv2+rural1+emp+sinsch+s1q5y+I(s1q5y^2)+female+hprim+hmid+hpmid+hwrk+rlh+hhsize+mih+fih+chr+mos+trd+as.factor(EcZone)+as.factor(region)"))
Mod.1 = regSARAR(formula = fml.1) #use full sample
summary.g2sls(Mod.1,marg_effect=FALSE,bet.id=3,L=1000,cl=5)
cat(rep("=",45)); cat("\n")
#----------------------------------------------------------------------------------
# Exclude religious dummies
cat("\n")
cat("Full Sample - Results with religious dummies excluded","\n")
fml.2 = as.formula(paste("wrk~pv2+rural1+emp+sinsch+s1q5y+I(s1q5y^2)+female+hprim+hmid+hpmid+hwrk+rlh+hhsize+mih+fih+as.factor(EcZone)+as.factor(region)"))
Mod.2 = regSARAR(formula = fml.2) #use full sample
summary.g2sls(Mod.2,marg_effect=FALSE,bet.id=3,L=1000,cl=5)
cat(rep("=",45)); cat("\n")
#----------------------------------------------------------------------------------
# Exclude household characteristics
cat("\n")
cat("Full Sample - Results with household characteristics excluded","\n")
fml.3 = as.formula(paste("wrk~pv2+rural1+emp+sinsch+s1q5y+I(s1q5y^2)+female+chr+mos+trd+as.factor(EcZone)+as.factor(region)"))
Mod.3 = regSARAR(formula = fml.3) #use full sample
summary.g2sls(Mod.3,marg_effect=FALSE,bet.id=3,L=1000,cl=5)
cat(rep("=",45)); cat("\n")
#----------------------------------------------------------------------------------
# Exclude Ecological Zone dummies
cat("\n")
cat("Full Sample - Results with Ecological Zone dummies excluded","\n")
fml.4 = as.formula(paste("wrk~pv2+rural1+emp+sinsch+s1q5y+I(s1q5y^2)+female+hprim+hmid+hpmid+hwrk+rlh+hhsize+mih+fih+chr+mos+trd+as.factor(region)"))
Mod.4 = regSARAR(formula = fml.4) #use full sample
summary.g2sls(Mod.4,marg_effect=FALSE,bet.id=3,L=1000,cl=5)
cat(rep("=",45)); cat("\n")
#----------------------------------------------------------------------------------
# Exclude Regional dummies
cat("\n")
cat("Full Sample - Results with Regional dummies excluded","\n")
fml.5 = as.formula(paste("wrk~pv2+rural1+emp+sinsch+s1q5y+I(s1q5y^2)+female+hprim+hmid+hpmid+hwrk+rlh+hhsize+mih+fih+chr+mos+trd+as.factor(EcZone)"))
Mod.5 = regSARAR(formula = fml.5) #use full sample
summary.g2sls(Mod.5,marg_effect=FALSE,bet.id=3,L=1000,cl=5)
cat(rep("=",45)); cat("\n")
#----------------------------------------------------------------------------------
# Exclude age and age squared
cat("\n")
cat("Full Sample - Results with age variables excluded","\n")
fml.8 = as.formula(paste("wrk~pv2+rural1+emp+sinsch+female+hprim+hmid+hpmid+hwrk+rlh+hhsize+mih+fih+chr+mos+trd+as.factor(EcZone)+as.factor(region)"))
Mod.8 = regSARAR(formula = fml.8) #use full sample
summary.g2sls(Mod.8,marg_effect=FALSE,bet.id=3,L=1000,cl=5)
cat(rep("=",45)); cat("\n")
#----------------------------------------------------------------------------------



#----------------------------------------------------------------------------------
# Formulae for sub-groups to study - Ecological Zones, Sex, Rural/Urban, Regions
fml.EcZone = as.formula(paste("wrk~pv2+rural1+emp+sinsch+s1q5y+I(s1q5y^2)+female+hprim+hmid+hpmid+hwrk+rlh+hhsize+mih+fih+chr+mos+trd+as.factor(region)"))
fml.sex = as.formula(paste("wrk~pv2+rural1+emp+sinsch+s1q5y+I(s1q5y^2)+hprim+hmid+hpmid+hwrk+rlh+hhsize+mih+fih+chr+mos+trd+as.factor(EcZone)+as.factor(region)"))
fml.rural = as.formula(paste("wrk~pv2+emp+sinsch+s1q5y+I(s1q5y^2)+female+hprim+hmid+hpmid+hwrk+rlh+hhsize+mih+fih+chr+mos+trd+as.factor(EcZone)+as.factor(region)"))
#for regions with more than one ecological zone
fml.region1 = as.formula(paste("wrk~pv2+rural1+emp+sinsch+s1q5y+I(s1q5y^2)+female+hprim+hmid+hpmid+hwrk+rlh+hhsize+mih+fih+chr+mos+trd+as.factor(EcZone)"))
#drop ecological zone dummies for regions with one ecological zone
fml.region2 = as.formula(paste("wrk~pv2+rural1+emp+sinsch+s1q5y+I(s1q5y^2)+female+hprim+hmid+hpmid+hwrk+rlh+hhsize+mih+fih+chr+mos+trd"))


# ------------- Ecological Zones ------------- 
cat("\n")
cat("Sub Sample Results by Ecological Zone","\n")
ssregs(var.sub=datk$EcZone,formula=fml.EcZone)

# ------------- Sex ------------- 
cat("\n")
cat("Sub Sample Results by Sex","\n")
ssregs(var.sub=datk$sex,formula=fml.sex)

# ------------- Rural/Urban ------------- 
cat("\n")
cat("Sub Sample Results by Rural/Urban (1/0)","\n")
ssregs(var.sub=datk$rural1,formula=fml.rural)
#----------------------------------------------------------------------------------



#----------------------------------------------------------------------------------
# ------------- Regions ------------- 
cat("\n")
cat("Sub Sample Results by Region","\n")
u.varR=unique(datk$region)

for(r in 1:length(unique(datk$region))){
  cat("Region Name: ",u.varR[r],"\n")
  cat("Ecological Zones in the ",u.varR[r]," Region are ")
  u.EZ=unique(datk$EcZone[which(datk$region==u.varR[r])])
  print(u.EZ); cat("\n")
  id = which(datk$region==u.varR[r])
  # avoid Ecological Zone dummies if there is only one Ecological in Region
  if(length(u.EZ)>1){
  Modl = regSARAR(formula = fml.region1,subsamp = id) 
  }else{
  Modl = regSARAR(formula = fml.region2,subsamp = id)
  }#end if/else
  print(summary.g2sls(Modl,marg_effect=FALSE,bet.id=3,L=1000,cl=5))
  cat(rep("=",45)); cat("\n")
}
#----------------------------------------------------------------------------------
# ------------- Region & Ecological Zone ------------- 
#u.eZ=unique(datk$EcZone)
for(r in 1:length(u.varR)){
  u.eZ=unique(datk$EcZone[which(datk$region==u.varR[r])])
  if(length(u.eZ)>1){#avoids redundancy where there is only one Ecological Zone in the Region
  for(e in 1:length(u.eZ)){
    cat("Region Name: ",u.varR[r]," & Ecological Zone: ",u.eZ[e],"\n")

    id = which(datk$region==u.varR[r] & datk$EcZone==u.eZ[e])
    Modl = regSARAR(formula = fml.region2,subsamp = id)
    print(summary.g2sls(Modl,marg_effect=FALSE,bet.id=3,L=1000,cl=5))
    cat(rep("=",45)); cat("\n")
  }#end for e
  }#end if
}#end for r
#----------------------------------------------------------------------------------

# List of Regions and their Ecological Zones
u.varR=unique(datk$region)

for(r in 1:length(u.varR)){
  cat("Region Name: ",u.varR[r],"\n")
  u.eZ=unique(datk$EcZone[which(datk$region==u.varR[r])])
  cat("Ecological Zones: ",u.eZ,"\n")
}


# Region Name:  Western 
# Ecological Zones:  Coastal Forest 
# Region Name:  Central 
# Ecological Zones:  Coastal Forest 
# Region Name:  Greater Accra 
# Ecological Zones:  Coastal 
# Region Name:  Volta 
# Ecological Zones:  Coastal Forest Savannah 
# Region Name:  Eastern 
# Ecological Zones:  Forest 
# Region Name:  Ashanti 
# Ecological Zones:  Forest 
# Region Name:  Brong Ahafo 
# Ecological Zones:  Savannah Forest 
# Region Name:  Northern 
# Ecological Zones:  Savannah 
# Region Name:  Upper East 
# Ecological Zones:  Savannah 
# Region Name:  Upper West 
# Ecological Zones:  Savannah

#=====================================================================================>
# Closing Log File
closeAllConnections() # Close connection to log file
#=====================================================================================>

Full Sample - Results 
Printing heteroskedasticity-consistent results...  
                                   Estimates Std. Error HC   t value HC  Pr(>|t|) HC
(Intercept)                    -0.5301589420  0.0397552530 -13.33556955 0.000000e+00
rho                             0.5478317046  0.0406081907  13.49067012 0.000000e+00
pv2                             0.0155772860  0.0057277107   2.71963561 6.535389e-03
rural1                          0.0303279288  0.0068290565   4.44101301 8.953638e-06
emp                             0.2834980117  0.0359312442   7.89001378 3.108624e-15
sinsch                         -0.1133523670  0.0224494517  -5.04922653 4.436025e-07
s1q5y                           0.0592223018  0.0041655953  14.21700810 0.000000e+00
I(s1q5y^2)                     -0.0014658966  0.0001953163  -7.50524292 6.128431e-14
female                         -0.0126025307  0.0046142928  -2.73119444 6.310523e-03
hprim                          -0.0358057142  0.0095120056  -3.76426546 1.670394e-04
hmid                           -0.0361870856  0.0063987802  -5.65531003 1.555652e-08
hpmid                          -0.0643558229  0.0078250922  -8.22428941 2.220446e-16
hwrk                            0.0606209370  0.0071170430   8.51771404 0.000000e+00
rlh                            -0.0023625681  0.0097086609  -0.24334644 8.077370e-01
hhsize                          0.0021341271  0.0008543066   2.49808096 1.248677e-02
mih                            -0.0069976947  0.0082509524  -0.84810751 3.963781e-01
fih                             0.0046539353  0.0069257067   0.67197984 5.015965e-01
chr                            -0.0032797846  0.0147638855  -0.22214915 8.241978e-01
mos                            -0.0131145571  0.0153759987  -0.85292392 3.937015e-01
trd                             0.0506275054  0.0258112783   1.96144898 4.982667e-02
as.factor(EcZone)Forest        -0.0004050522  0.0080443198  -0.05035258 9.598414e-01
as.factor(EcZone)Savannah       0.0006404662  0.0140326582   0.04564112 9.635963e-01
as.factor(region)Brong Ahafo    0.0396075296  0.0142578495   2.77794556 5.470378e-03
as.factor(region)Central       -0.0289937862  0.0104988649  -2.76161152 5.751687e-03
as.factor(region)Eastern        0.0322226447  0.0112808822   2.85639404 4.284830e-03
as.factor(region)Greater Accra -0.0180948261  0.0122763574  -1.47395726 1.404931e-01
as.factor(region)Northern      -0.0096230732  0.0164165847  -0.58617997 5.577546e-01
as.factor(region)Upper East     0.0848862973  0.0179266811   4.73519313 2.188466e-06
as.factor(region)Upper West     0.0514393329  0.0171437081   3.00047882 2.695555e-03
as.factor(region)Volta         -0.0144507635  0.0114497232  -1.26210592 2.069107e-01
as.factor(region)Western       -0.0062186024  0.0110966993  -0.56040109 5.752059e-01

 For goodness-of-fit measures viz. R-squared, Wald Tests. 

Call:
AER::ivreg(formula = y ~ as.matrix(R[, -1]) | as.matrix(Z[, -1]), 
    weights = 1/e2)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.7300 -0.4330 -0.1873  0.5894  2.5323 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.53016    0.01714  -30.92   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.5325 on 21174 degrees of freedom
Multiple R-Squared: 0.7363,	Adjusted R-squared: 0.7359 
Wald test:  1361 on 30 and 21174 DF,  p-value: < 2.2e-16 



 Lambda =  0.4946998 
= = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =

Full Sample - Results with religious dummies excluded 
Printing heteroskedasticity-consistent results...  
                                   Estimates Std. Error HC   t value HC  Pr(>|t|) HC
(Intercept)                    -0.5483987343  0.0371932776 -14.74456595 0.000000e+00
rho                             0.5171152502  0.0410728840  12.59018604 0.000000e+00
pv2                             0.0177928673  0.0057290589   3.10572254 1.898147e-03
rural1                          0.0312870873  0.0068854359   4.54395156 5.520935e-06
emp                             0.2911038694  0.0362130220   8.03865166 8.881784e-16
sinsch                         -0.1178209378  0.0224612915  -5.24551039 1.558503e-07
s1q5y                           0.0581335320  0.0041671891  13.95029859 0.000000e+00
I(s1q5y^2)                     -0.0013603495  0.0001953594  -6.96331811 3.323564e-12
female                         -0.0131854364  0.0046149028  -2.85714284 4.274734e-03
hprim                          -0.0349446152  0.0094982262  -3.67906747 2.340883e-04
hmid                           -0.0274774157  0.0063044800  -4.35839528 1.310196e-05
hpmid                          -0.0617068788  0.0076513196  -8.06486748 6.661338e-16
hwrk                            0.0642337981  0.0071218727   9.01922866 0.000000e+00
rlh                            -0.0036866055  0.0097018004  -0.37999189 7.039514e-01
hhsize                          0.0020309206  0.0008506356   2.38753301 1.696188e-02
mih                            -0.0054048999  0.0082375488  -0.65612964 5.117407e-01
fih                            -0.0005063043  0.0069095625  -0.07327588 9.415866e-01
as.factor(EcZone)Forest         0.0029024305  0.0080429634   0.36086581 7.181998e-01
as.factor(EcZone)Savannah       0.0016692253  0.0140227162   0.11903723 9.052459e-01
as.factor(region)Brong Ahafo    0.0494852699  0.0142870366   3.46364828 5.329028e-04
as.factor(region)Central       -0.0270168706  0.0105099600  -2.57059690 1.015234e-02
as.factor(region)Eastern        0.0350050953  0.0113236427   3.09132815 1.992633e-03
as.factor(region)Greater Accra -0.0087691986  0.0122820391  -0.71398556 4.752361e-01
as.factor(region)Northern       0.0005275896  0.0162315599   0.03250394 9.740702e-01
as.factor(region)Upper East     0.1063353303  0.0178562700   5.95506959 2.599612e-09
as.factor(region)Upper West     0.0639129085  0.0170463687   3.74935622 1.772891e-04
as.factor(region)Volta         -0.0023587678  0.0114536078  -0.20594103 8.368370e-01
as.factor(region)Western        0.0038348190  0.0111007957   0.34545443 7.297528e-01

 For goodness-of-fit measures viz. R-squared, Wald Tests. 

Call:
AER::ivreg(formula = y ~ as.matrix(R[, -1]) | as.matrix(Z[, -1]), 
    weights = 1/e2)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.2648 -0.4368 -0.1861  0.6144  2.1170 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.54840    0.01654  -33.15   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.5402 on 21177 degrees of freedom
Multiple R-Squared: 0.7349,	Adjusted R-squared: 0.7346 
Wald test:  1499 on 27 and 21177 DF,  p-value: < 2.2e-16 



 Lambda =  0.5623985 
= = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =

Full Sample - Results with household characteristics excluded 
Printing heteroskedasticity-consistent results...  
                                   Estimates Std. Error HC  t value HC  Pr(>|t|) HC
(Intercept)                    -0.5754485011   0.044080932 -13.0543632 0.000000e+00
rho                             0.0991254148   0.059844582   1.6563808 9.764473e-02
pv2                             0.0022303427   0.005794987   0.3848745 7.003304e-01
rural1                          0.0637303593   0.008234688   7.7392564 9.992007e-15
emp                             0.5325188848   0.049702725  10.7140783 0.000000e+00
sinsch                         -0.1351352766   0.022766308  -5.9357572 2.924925e-09
s1q5y                           0.0468926510   0.004180649  11.2165973 0.000000e+00
I(s1q5y^2)                     -0.0004070956   0.000196249  -2.0743830 3.804375e-02
female                         -0.0149126137   0.004634456  -3.2177699 1.291914e-03
chr                            -0.0035343619   0.014713179  -0.2402174 8.101617e-01
mos                            -0.0256582233   0.015397864  -1.6663495 9.564383e-02
trd                             0.1088561903   0.031000351   3.5114503 4.456688e-04
as.factor(EcZone)Forest        -0.0207710966   0.008096022  -2.5655929 1.029996e-02
as.factor(EcZone)Savannah       0.0072870082   0.014224872   0.5122723 6.084604e-01
as.factor(region)Brong Ahafo    0.0794920817   0.015117902   5.2581424 1.455178e-07
as.factor(region)Central       -0.0485612342   0.011045257  -4.3965688 1.099755e-05
as.factor(region)Eastern        0.0290104930   0.011873465   2.4433047 1.455344e-02
as.factor(region)Greater Accra -0.0701123563   0.012572055  -5.5768415 2.449249e-08
as.factor(region)Northern      -0.0189284993   0.016502938  -1.1469775 2.513909e-01
as.factor(region)Upper East     0.1513155818   0.019062193   7.9379943 1.998401e-15
as.factor(region)Upper West     0.0299047212   0.017895307   1.6710930 9.470331e-02
as.factor(region)Volta         -0.0328636309   0.011516533  -2.8536046 4.322630e-03
as.factor(region)Western       -0.0151338955   0.011160721  -1.3559962 1.751004e-01

 For goodness-of-fit measures viz. R-squared, Wald Tests. 

Call:
AER::ivreg(formula = y ~ as.matrix(R[, -1]) | as.matrix(Z[, -1]), 
    weights = 1/e2)

Residuals:
    Min      1Q  Median      3Q     Max 
-5.5086 -0.4945 -0.1853  0.7461 10.7104 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.57545    0.02207  -26.07   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.7009 on 21182 degrees of freedom
Multiple R-Squared: 0.5276,	Adjusted R-squared: 0.5271 
Wald test: 933.4 on 22 and 21182 DF,  p-value: < 2.2e-16 



 Lambda =  0.8770424 
= = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =

Full Sample - Results with Ecological Zone dummies excluded 
Printing heteroskedasticity-consistent results...  
                                   Estimates Std. Error HC   t value HC  Pr(>|t|) HC
(Intercept)                    -0.5574054714  0.0391481836 -14.23834826 0.000000e+00
rho                             0.5721088487  0.0407839495  14.02779415 0.000000e+00
pv2                             0.0169943840  0.0057334739   2.96406407 3.036051e-03
rural1                          0.0320995039  0.0068337143   4.69722652 2.637179e-06
emp                             0.2693621555  0.0361408015   7.45313176 9.126033e-14
sinsch                         -0.1074864262  0.0224428376  -4.78934207 1.673290e-06
s1q5y                           0.0642028940  0.0041653482  15.41357174 0.000000e+00
I(s1q5y^2)                     -0.0016911688  0.0001953165  -8.65860658 0.000000e+00
female                         -0.0116381833  0.0046141126  -2.52230152 1.165897e-02
hprim                          -0.0360858404  0.0094865157  -3.80390880 1.424306e-04
hmid                           -0.0275402350  0.0063773912  -4.31841708 1.571522e-05
hpmid                          -0.0544977593  0.0078065836  -6.98099993 2.930767e-12
hwrk                            0.0595265410  0.0071113792   8.37060420 0.000000e+00
rlh                            -0.0030209466  0.0097056224  -0.31125738 7.556050e-01
hhsize                          0.0016645002  0.0008539474   1.94918356 5.127351e-02
mih                            -0.0103344986  0.0082506212  -1.25257219 2.103614e-01
fih                             0.0032513762  0.0069252721   0.46949436 6.387163e-01
chr                            -0.0040490706  0.0147675750  -0.27418656 7.839413e-01
mos                            -0.0070918228  0.0153739355  -0.46128871 6.445915e-01
trd                             0.0462700816  0.0249924471   1.85136259 6.411741e-02
as.factor(region)Brong Ahafo    0.0438908175  0.0119604739   3.66965540 2.428777e-04
as.factor(region)Central       -0.0187694275  0.0098341079  -1.90860499 5.631307e-02
as.factor(region)Eastern        0.0362463484  0.0112821172   3.21272574 1.314818e-03
as.factor(region)Greater Accra -0.0105641374  0.0095597758  -1.10506120 2.691331e-01
as.factor(region)Northern      -0.0001032071  0.0111561047  -0.00925118 9.926187e-01
as.factor(region)Upper East     0.0932686612  0.0132412566   7.04379231 1.870726e-12
as.factor(region)Upper West     0.0586144702  0.0122186861   4.79711726 1.609654e-06
as.factor(region)Volta         -0.0060518099  0.0105673731  -0.57268820 5.668558e-01
as.factor(region)Western        0.0051848502  0.0103294303   0.50194929 6.157032e-01

 For goodness-of-fit measures viz. R-squared, Wald Tests. 

Call:
AER::ivreg(formula = y ~ as.matrix(R[, -1]) | as.matrix(Z[, -1]), 
    weights = 1/e2)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.6789 -0.4347 -0.1960  0.5606  2.0394 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.55741    0.01634  -34.11   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.5285 on 21176 degrees of freedom
Multiple R-Squared: 0.7631,	Adjusted R-squared: 0.7628 
Wald test:  1737 on 28 and 21176 DF,  p-value: < 2.2e-16 



 Lambda =  0.4326917 
= = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =

Full Sample - Results with Regional dummies excluded 
Printing heteroskedasticity-consistent results...  
                             Estimates Std. Error HC  t value HC  Pr(>|t|) HC
(Intercept)               -0.531795137  0.0384954271 -13.8145015 0.000000e+00
rho                        0.770967292  0.0365024448  21.1209769 0.000000e+00
pv2                        0.014052269  0.0057078525   2.4619188 1.381960e-02
rural1                     0.014608367  0.0065082554   2.2445903 2.479445e-02
emp                        0.171690247  0.0326726699   5.2548582 1.481386e-07
sinsch                    -0.074797438  0.0226781539  -3.2982155 9.730144e-04
s1q5y                      0.069320550  0.0041609739  16.6596937 0.000000e+00
I(s1q5y^2)                -0.002151049  0.0001953767 -11.0097513 0.000000e+00
female                    -0.004259037  0.0046126540  -0.9233377 3.558312e-01
hprim                     -0.029340243  0.0094837126  -3.0937508 1.976435e-03
hmid                      -0.030624791  0.0063134318  -4.8507359 1.230042e-06
hpmid                     -0.057656696  0.0077031602  -7.4848107 7.172041e-14
hwrk                       0.058312343  0.0070982003   8.2150885 2.220446e-16
rlh                       -0.009601264  0.0097021624  -0.9896004 3.223695e-01
hhsize                     0.001435886  0.0008380845   1.7132952 8.665825e-02
mih                       -0.003832785  0.0082384840  -0.4652294 6.417672e-01
fih                        0.009626524  0.0069064571   1.3938440 1.633647e-01
chr                       -0.005123800  0.0147320523  -0.3477995 7.279908e-01
mos                       -0.018627512  0.0153347253  -1.2147275 2.244701e-01
trd                        0.020587482  0.0236543291   0.8703473 3.841106e-01
as.factor(EcZone)Forest    0.011688430  0.0063985328   1.8267360 6.773947e-02
as.factor(EcZone)Savannah  0.025849666  0.0088610847   2.9172124 3.531752e-03

 For goodness-of-fit measures viz. R-squared, Wald Tests. 

Call:
AER::ivreg(formula = y ~ as.matrix(R[, -1]) | as.matrix(Z[, -1]), 
    weights = 1/e2)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.1572 -0.4051 -0.1571  0.4227  2.6311 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.53180    0.01506  -35.31   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.509 on 21183 degrees of freedom
Multiple R-Squared: 0.7535,	Adjusted R-squared: 0.7533 
Wald test:  1908 on 21 and 21183 DF,  p-value: < 2.2e-16 



 Lambda =  -0.6386515 
= = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =

Full Sample - Results with age variables excluded 
Printing heteroskedasticity-consistent results...  
                                   Estimates Std. Error HC  t value HC  Pr(>|t|) HC
(Intercept)                     0.0723354242  0.0374498417  1.93152817 5.341777e-02
rho                             0.8637856962  0.0498975023 17.31120107 0.000000e+00
pv2                             0.0053310522  0.0061191670  0.87120555 3.836419e-01
rural1                          0.0044711453  0.0075771450  0.59008311 5.551349e-01
emp                             0.0681859529  0.0429497516  1.58757502 1.123825e-01
sinsch                         -0.1336819695  0.0228319370 -5.85504285 4.768870e-09
female                         -0.0027352832  0.0048785750 -0.56067257 5.750208e-01
hprim                          -0.0205591691  0.0100829948 -2.03899432 4.145059e-02
hmid                           -0.0153447104  0.0068652019 -2.23514335 2.540794e-02
hpmid                          -0.0230172084  0.0081878754 -2.81113320 4.936734e-03
hwrk                            0.0384351243  0.0073521989  5.22770465 1.716274e-07
rlh                             0.0195897864  0.0102455755  1.91202402 5.587311e-02
hhsize                          0.0012129231  0.0009233501  1.31361129 1.889770e-01
mih                            -0.0219459557  0.0086410725 -2.53972593 1.109394e-02
fih                            -0.0111787517  0.0073062052 -1.53003527 1.260080e-01
chr                            -0.0046582486  0.0159160602 -0.29267598 7.697698e-01
mos                            -0.0069620737  0.0165651112 -0.42028536 6.742770e-01
trd                             0.0076488115  0.0848034296  0.09019460 9.281326e-01
as.factor(EcZone)Forest        -0.0037987957  0.0082600173 -0.45990166 6.455868e-01
as.factor(EcZone)Savannah      -0.0065199857  0.0151433739 -0.43055040 6.667953e-01
as.factor(region)Brong Ahafo    0.0074747724  0.0155928406  0.47937208 6.316740e-01
as.factor(region)Central       -0.0014792898  0.0110301925 -0.13411278 8.933134e-01
as.factor(region)Eastern        0.0075720066  0.0121978808  0.62076411 5.347548e-01
as.factor(region)Greater Accra  0.0006122422  0.0125258778  0.04887818 9.610164e-01
as.factor(region)Northern      -0.0104805401  0.0175985830 -0.59553318 5.514871e-01
as.factor(region)Upper East     0.0339120779  0.0196591114  1.72500563 8.452646e-02
as.factor(region)Upper West     0.0123781921  0.0186295364  0.66443908 5.064093e-01
as.factor(region)Volta         -0.0012529309  0.0119808783 -0.10457755 9.167110e-01
as.factor(region)Western       -0.0062602934  0.0117422335 -0.53314332 5.939344e-01

 For goodness-of-fit measures viz. R-squared, Wald Tests. 

Call:
AER::ivreg(formula = y ~ as.matrix(R[, -1]) | as.matrix(Z[, -1]), 
    weights = 1/e2)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.29168 -0.38192 -0.04063  0.34819  2.13240 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.07234    0.01469   4.923 8.57e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.4958 on 21176 degrees of freedom
Multiple R-Squared: 0.6388,	Adjusted R-squared: 0.6384 
Wald test: 656.2 on 28 and 21176 DF,  p-value: < 2.2e-16 



 Lambda =  -0.999959 
= = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =

Sub Sample Results by Ecological Zone 
Sub-sample is  Coastal  
Error in base::chol2inv(x, ...) : 
  element (52, 52) is zero, so the inverse cannot be computed
Called from: base::chol2inv(x, ...)

Sub Sample Results by Sex 
Sub-sample is  Female  
Error in base::chol2inv(x, ...) : 
  element (57, 57) is zero, so the inverse cannot be computed
Called from: base::chol2inv(x, ...)

Sub Sample Results by Rural/Urban (1/0) 
Sub-sample is  0  
Printing heteroskedasticity-consistent results...  
                                   Estimates Std. Error HC t value HC  Pr(>|t|) HC
(Intercept)                    -0.1608970750   0.060468288 -2.6608505 0.0077943554
rho                             0.1756360465   0.053291857  3.2957389 0.0009816320
pv2                            -0.0037632735   0.010812709 -0.3480417 0.7278088661
emp                             0.1270945440   0.037669700  3.3739197 0.0007410597
sinsch                         -0.0591835873   0.032570068 -1.8171159 0.0691993715
s1q5y                           0.0152116319   0.005590957  2.7207563 0.0065132763
I(s1q5y^2)                     -0.0003624094   0.000265496 -1.3650276 0.1722443875
female                          0.0069148230   0.006478165  1.0674046 0.2857891587
hprim                          -0.0158309818   0.014335651 -1.1043085 0.2694593212
hmid                           -0.0190963825   0.009492182 -2.0118011 0.0442409129
hpmid                          -0.0358536073   0.009720366 -3.6885036 0.0002255768
hwrk                            0.0192352718   0.008581346  2.2415216 0.0249923109
rlh                             0.0065321637   0.013115905  0.4980338 0.6184602353
hhsize                          0.0043778940   0.001573332  2.7825617 0.0053931604
mih                            -0.0088352900   0.010710814 -0.8248944 0.4094315646
fih                            -0.0048928314   0.009248334 -0.5290501 0.5967707289
chr                             0.0046994900   0.031838053  0.1476061 0.8826536602
mos                            -0.0059787657   0.033086307 -0.1807021 0.8566013975
trd                            -0.0075816246   0.037527538 -0.2020283 0.8398946008
as.factor(EcZone)Forest         0.0026831820   0.010501374  0.2555077 0.7983309930
as.factor(EcZone)Savannah       0.0026210760   0.022666564  0.1156362 0.9079408471
as.factor(region)Brong Ahafo    0.0266395113   0.020134157  1.3231004 0.1858019969
as.factor(region)Central       -0.0040164481   0.012369381 -0.3247089 0.7454014018
as.factor(region)Eastern        0.0335613850   0.015706536  2.1367782 0.0326160323
as.factor(region)Greater Accra -0.0077536996   0.014449145 -0.5366200 0.5915301453
as.factor(region)Northern       0.0084195415   0.026207586  0.3212635 0.7480106922
as.factor(region)Upper East     0.0581075003   0.030322477  1.9163177 0.0553246641
as.factor(region)Upper West    -0.0096594337   0.026019647 -0.3712361 0.7104616505
as.factor(region)Volta         -0.0070601645   0.014667688 -0.4813413 0.6302739215
as.factor(region)Western       -0.0045832133   0.014518559 -0.3156796 0.7522456878

 For goodness-of-fit measures viz. R-squared, Wald Tests. 

Call:
AER::ivreg(formula = y ~ as.matrix(R[, -1]) | as.matrix(Z[, -1]), 
    weights = 1/e2)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.42015 -0.16609 -0.09441  0.11227  1.78059 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.16090    0.02064  -7.796 7.19e-15 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.4419 on 7972 degrees of freedom
Multiple R-Squared: 0.2156,	Adjusted R-squared: 0.2127 
Wald test: 30.17 on 29 and 7972 DF,  p-value: < 2.2e-16 



 Lambda =  0.8454247 
= = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
Sub-sample is  1  
Printing heteroskedasticity-consistent results...  
                                   Estimates Std. Error HC   t value HC  Pr(>|t|) HC
(Intercept)                    -0.7376764455  0.0505190426 -14.60194825 0.000000e+00
rho                             0.5933677465  0.0460983800  12.87177004 0.000000e+00
pv2                             0.0119024225  0.0067190293   1.77144970 7.648595e-02
emp                             0.2956603175  0.0455189464   6.49532427 8.285483e-11
sinsch                         -0.0913871637  0.0304394423  -3.00226144 2.679819e-03
s1q5y                           0.0777710609  0.0056827674  13.68542046 0.000000e+00
I(s1q5y^2)                     -0.0016885330  0.0002661022  -6.34542997 2.218050e-10
female                         -0.0178483610  0.0062075505  -2.87526634 4.036868e-03
hprim                          -0.0305526158  0.0123455838  -2.47478096 1.333179e-02
hmid                           -0.0315760755  0.0084627561  -3.73118111 1.905841e-04
hpmid                          -0.0467872250  0.0136753219  -3.42128874 6.232512e-04
hwrk                            0.0968052006  0.0108432297   8.92770908 0.000000e+00
rlh                            -0.0081951185  0.0132638413  -0.61785409 5.366715e-01
hhsize                         -0.0002271539  0.0009949606  -0.22830442 8.194096e-01
mih                            -0.0091054835  0.0114574438  -0.79472208 4.267752e-01
fih                             0.0158082414  0.0099104644   1.59510601 1.106885e-01
chr                            -0.0066743659  0.0164676861  -0.40530077 6.852564e-01
mos                            -0.0086078180  0.0173037986  -0.49745251 6.188700e-01
trd                             0.0664292399  0.0480611345   1.38218210 1.669158e-01
as.factor(EcZone)Forest         0.0011256741  0.0114850318   0.09801227 9.219225e-01
as.factor(EcZone)Savannah      -0.0010750797  0.0177252879  -0.06065232 9.516361e-01
as.factor(region)Brong Ahafo    0.0749701530  0.0197980743   3.78673965 1.526369e-04
as.factor(region)Central       -0.0270811773  0.0160569159  -1.68657403 9.168529e-02
as.factor(region)Eastern        0.0408398322  0.0157343796   2.59557944 9.443157e-03
as.factor(region)Greater Accra  0.0306485749  0.0282639761   1.08436884 2.782013e-01
as.factor(region)Northern       0.0083711738  0.0217351144   0.38514514 7.001299e-01
as.factor(region)Upper East     0.1035183918  0.0231018016   4.48096619 7.430588e-06
as.factor(region)Upper West     0.0928364519  0.0226079080   4.10637074 4.019242e-05
as.factor(region)Volta          0.0056350810  0.0166191492   0.33907157 7.345558e-01
as.factor(region)Western        0.0235103579  0.0159778067   1.47143837 1.411726e-01

 For goodness-of-fit measures viz. R-squared, Wald Tests. 

Call:
AER::ivreg(formula = y ~ as.matrix(R[, -1]) | as.matrix(Z[, -1]), 
    weights = 1/e2)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.0530 -0.5047 -0.1197  0.5906  1.8242 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.73768    0.02121  -34.78   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.5651 on 13173 degrees of freedom
Multiple R-Squared: 0.9537,	Adjusted R-squared: 0.9536 
Wald test:  8866 on 29 and 13173 DF,  p-value: < 2.2e-16 



 Lambda =  0.4601361 
= = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =

Sub Sample Results by Region 
Region Name:  Western 
Ecological Zones in the  Western  Region are [1] "Coastal" "Forest" 

Error in base::chol2inv(x, ...) : 
  element (50, 50) is zero, so the inverse cannot be computed
Called from: base::chol2inv(x, ...)
debug at #2: u.eZ = unique(datk$EcZone[which(datk$region == u.varR[r])])
debug at #3: if (length(u.eZ) > 1) {
    for (e in 1:length(u.eZ)) {
        cat("Region Name: ", u.varR[r], " & Ecological Zone: ", 
            u.eZ[e], "\n")
        id = which(datk$region == u.varR[r] & datk$EcZone == 
            u.eZ[e])
        Modl = regSARAR(formula = fml.region2, subsamp = id)
        print(summary.g2sls(Modl, marg_effect = FALSE, bet.id = 3, 
            L = 1000, cl = 5))
        cat(rep("=", 45))
        cat("\n")
    }
}
debug at #4: for (e in 1:length(u.eZ)) {
    cat("Region Name: ", u.varR[r], " & Ecological Zone: ", u.eZ[e], 
        "\n")
    id = which(datk$region == u.varR[r] & datk$EcZone == u.eZ[e])
    Modl = regSARAR(formula = fml.region2, subsamp = id)
    print(summary.g2sls(Modl, marg_effect = FALSE, bet.id = 3, 
        L = 1000, cl = 5))
    cat(rep("=", 45))
    cat("\n")
}
debug at #2: cat("Region Name: ", u.varR[r], "\n")
Region Name:  Western 
debug at #3: u.eZ = unique(datk$EcZone[which(datk$region == u.varR[r])])
debug at #4: cat("Ecological Zones: ", u.eZ, "\n")
Ecological Zones:  Coastal Forest 
debug at #2: cat("Region Name: ", u.varR[r], "\n")
