Does early planting of wheat affect mean yield and risk? Just-Pope residual based and Antle moments production risk analysis

1 Introduction

There is ample evidence that early sowing of long duration wheat varieties increases yieCSISA_KVKestim and profitability of wheat production in EIGP. However, there are concerns that it may affect the riskiness of the production system as well as alter the effectiveness of other inputs like fertilizers and irrigation. We use the residual based and moments-based approaches to econometric assessment of productivity and risk of different sowing date strategies and variety maturity class. We also investigate the complementary and substitution patterns between these inputs and other production inputs including weather.

This notebook provides a reproducible workflow for Just-Pope production function and moments based approach to production risk analysis.

We use data collected under the CSISA-KVK agronomic trials.

CSISA_KVKestim=read.csv("CSISA_KVK_wheat_public_cleaned.csv")

#load("CSISA_KVK_Public_Workspace.RData")

table(CSISA_KVKestim$SowingSchedule) 

 T1  T2  T3  T4  T5 
408 858 806 881 695 

#CSISA_KVKestim$SowingSchedule=ordered(CSISA_KVKestim$SowingSchedule,levels=c("T5","T4","T3","T2","T1"))

#Create dummies for some categorical variables
library(fastDummies)
CSISA_KVKestim <- fastDummies::dummy_cols(CSISA_KVKestim, select_columns =c("SowingSchedule","VarietyClass","SoilType","CropEstablishment","Year","District"))

# Make calculations on nitrogen, phosphurus, potassium applied 

## Nitrogen
CSISA_KVKestim$GradeNPKN[CSISA_KVKestim$GradeNPK=="10.26.26"]=0.10
CSISA_KVKestim$GradeNPKN[CSISA_KVKestim$GradeNPK %in% c("12.32.16", "12:32:16", "12:32:16 PM")] <- 0.12
CSISA_KVKestim$GradeNPKN[CSISA_KVKestim$GradeNPK == "14:35:14"] <- 0.14
CSISA_KVKestim$GradeNPKN[CSISA_KVKestim$GradeNPK == "18:46:00"] <- 0.18
CSISA_KVKestim$GradeNPKN[CSISA_KVKestim$GradeNPK == "20.20.0.13"] <- 0.20

## Phosphorus
CSISA_KVKestim$GradeNPKP[CSISA_KVKestim$GradeNPK == "10.26.26"] <- 0.26
CSISA_KVKestim$GradeNPKP[CSISA_KVKestim$GradeNPK %in% c("12.32.16", "12:32:16", "12:32:16 PM")] <- 0.32
CSISA_KVKestim$GradeNPKP[CSISA_KVKestim$GradeNPK == "14:35:14"] <- 0.35
CSISA_KVKestim$GradeNPKP[CSISA_KVKestim$GradeNPK == "18:46:00"] <- 0.46
CSISA_KVKestim$GradeNPKP[CSISA_KVKestim$GradeNPK == "20.20.0.13"] <- 0.20

# Potassium
CSISA_KVKestim$GradeNPKK[CSISA_KVKestim$GradeNPK == "10.26.26"] <- 0.26
CSISA_KVKestim$GradeNPKK[CSISA_KVKestim$GradeNPK %in% c("12.32.16", "12:32:16", "12:32:16 PM")] <- 0.16
CSISA_KVKestim$GradeNPKK[CSISA_KVKestim$GradeNPK == "14:35:14"] <- 0.14
CSISA_KVKestim$GradeNPKK[CSISA_KVKestim$GradeNPK == "18:46:00"] <- 0
CSISA_KVKestim$GradeNPKK[CSISA_KVKestim$GradeNPK == "20.20.0.13"] <- 0.13

# # 
# [21] "BasalDAP"
#  [22] "BasalNPK"
#  [23] "GradeNPK"
#  [24] "BasalMOP"
#  [25] "BasalZn"
#  [26] "Split1Urea"
#  [27] "Split2Urea"
#  [28] "Split3Urea"

# Nutrient Content ----------------------
# Taken from Cedrez, Chamberlain, Guo and Hijmans, p3
### N -----------------------------------
# CSISA_KVKestim$BasalDAPN=CSISA_KVKestim$BasalDAP*0.18 
# CSISA_KVKestim$Split1UreaN=CSISA_KVKestim$Split1Urea*0.46
# CSISA_KVKestim$Split2UreaN=CSISA_KVKestim$Split2Urea*0.46
# CSISA_KVKestim$Split3UreaN=CSISA_KVKestim$Split3Urea*0.46
# 
# CSISA_KVKestim$BasalNPKN=CSISA_KVKestim$BasalNPK*CSISA_KVKestim$GradeNPKN
# 
# CSISA_KVKestim$N=rowSums(CSISA_KVKestim[,c("BasalDAPN","Split1UreaN","Split2UreaN","Split3UreaN","BasalNPKN")],na.rm = TRUE)
# 
# CSISA_KVKestim$Nperha=CSISA_KVKestim$N/CSISA_KVKestim$C.q306_cropLarestAreaHA
# CSISA_KVKestim$NperhaSq=CSISA_KVKestim$Nperha*CSISA_KVKestim$Nperha
# 
# ### P ------------------------------------
# CSISA_KVKestim$F.totAmtDAPP=CSISA_KVKestim$F.totAmtDAP*0.46
# CSISA_KVKestim$F.totAmtUreaP=CSISA_KVKestim$F.totAmtUrea*0
# CSISA_KVKestim$F.totAmtNPKP=CSISA_KVKestim$F.totAmtNPK*CSISA_KVKestim$F.q51071_gradeNPKP
# 
# 
# CSISA_KVKestim$P2O5=rowSums(CSISA_KVKestim[,c("F.totAmtDAPP","F.totAmtUreaP","F.totAmtNPKP","F.totAmtTSPP","F.totAmtSSPP","F.totAmtNPKSP")],na.rm = TRUE)
# CSISA_KVKestim$P2O5perha=CSISA_KVKestim$P2O5/CSISA_KVKestim$C.q306_cropLarestAreaHA
# 
# 
# 
# # Get CHIRPS temperature and rainfall data(
# 
# # install.packages("remotes")
# # remotes::install_github("environmentalinformatics-marburg/heavyRain")
# 
# library("chirps")
# library("terra")
# 
# # Case 1: return as a data.frame
# dates <- c("2017-12-15","2017-12-31")
# lonlat <- data.frame(lon = c(-55.0281,-54.9857), lat = c(-2.8094, -2.8756))
# 
# r1 <- get_chirps(lonlat, dates, server = "CHC")

1.1 Graphics

# Bar graphs showing percentage of farmers adopting these practices  

library(tidyverse) 
library(ggplot2)  

bar_chart=function(dat,var){   dat|>     drop_na({{var}})|>     mutate({{var}}:=factor({{var}})|>fct_infreq())|>     ggplot()+     geom_bar(aes(y={{var}}),fill="dodgerblue4")+     theme_minimal(base_size = 16) }   

sow_plot=bar_chart(CSISA_KVKestim,SowingSchedule)+labs(y="Sowing dates") 
 
sow_plot

library(ggpubr) 
library(tidyverse) 

#Sowing dates 

SowingDate_Options_Errorplot=   CSISA_KVKestim%>%   
  drop_na(SowingSchedule) %>%   
  ggerrorplot(x = "SowingSchedule", y = "GrainYield",add = "mean", error.plot = "errorbar", color="steelblue", ggtheme=theme_bw())+   
  labs(x="Sowing date options",y="Wheat yield (t/ha)")+   
  theme_bw(base_size = 16)+
  coord_flip()  
SowingDate_Options_Errorplot 

1.2 Descriptives

library(fBasics)

summ_stats <- fBasics::basicStats(CSISA_KVKestim[,c("GrainYield","VarietyClass_LDV","SoilType_Heavy","SoilType_Medium","SoilType_Low","CropEstablishment_CT","CropEstablishment_CT-line","CropEstablishment_ZT")]) 

summ_stats <- as.data.frame(t(summ_stats)) 

# Rename some of the columns for convenience 

summ_stats <- summ_stats[c("Mean", "Stdev", "Minimum", "1. Quartile", "Median",  "3. Quartile", "Maximum")] %>%   
rename("Lower quartile" = '1. Quartile', "Upper quartile"= "3. Quartile")  

summ_stats 

                              Mean    Stdev Minimum Lower quartile Median
GrainYield                4.067177 1.035497    1.39           3.23   4.06
VarietyClass_LDV          0.710526 0.453580    0.00           0.00   1.00
SoilType_Heavy            0.214090 0.410246    0.00           0.00   0.00
SoilType_Medium           0.763980 0.424693    0.00           1.00   1.00
SoilType_Low              0.021930 0.146475    0.00           0.00   0.00
CropEstablishment_CT      0.020285 0.140993    0.00           0.00   0.00
CropEstablishment_CT.line 0.004386 0.066090    0.00           0.00   0.00
CropEstablishment_ZT      0.975329 0.155142    0.00           1.00   1.00
                          Upper quartile Maximum
GrainYield                          4.84    6.74
VarietyClass_LDV                    1.00    1.00
SoilType_Heavy                      0.00    1.00
SoilType_Medium                     1.00    1.00
SoilType_Low                        0.00    1.00
CropEstablishment_CT                0.00    1.00
CropEstablishment_CT.line           0.00    1.00
CropEstablishment_ZT                1.00    1.00

2 Just-pope production function

Just and Pope (1979) proposed the a three step estimation framework of the effect of inputs on mean yield and risk. The first step involves estimating a production function of any functional form (e.g., quadratic, cobb-douglas) using OLS. We then collect the residuals, square them and put them in a logarithm. We use this log (res^2) as dependent variable in the second stage estimation. This is the variance model. In the final step, we use the inverse of the exponential of the fitted values from the second stage as weights in a weighted least squares (WLS). This three step procedure is also called feasible generalized least sqaures (FGLS) estimation.

ols = lm(GrainYield ~ SowingSchedule+VarietyClass_LDV + SeedRate + CropEstablishment + IrrigationNumber + SoilType + LandType + as.factor(Year)+ District, data = CSISA_KVKestim)
summary(ols)

Call:
lm(formula = GrainYield ~ SowingSchedule + VarietyClass_LDV + 
    SeedRate + CropEstablishment + IrrigationNumber + SoilType + 
    LandType + as.factor(Year) + District, data = CSISA_KVKestim)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.56769 -0.38434 -0.00054  0.40224  2.43981 

Coefficients:
                          Estimate Std. Error t value Pr(>|t|)    
(Intercept)               2.772121   0.166600  16.639  < 2e-16 ***
SowingScheduleT2         -0.300546   0.038204  -7.867 4.76e-15 ***
SowingScheduleT3         -0.640101   0.039523 -16.196  < 2e-16 ***
SowingScheduleT4         -1.038753   0.041335 -25.130  < 2e-16 ***
SowingScheduleT5         -1.353139   0.043525 -31.089  < 2e-16 ***
VarietyClass_LDV          0.441173   0.025483  17.313  < 2e-16 ***
SeedRate                 -0.005214   0.002994  -1.741 0.081713 .  
CropEstablishmentCT-line  0.641219   0.173942   3.686 0.000231 ***
CropEstablishmentZT       0.662518   0.077162   8.586  < 2e-16 ***
IrrigationNumber          0.492991   0.018425  26.756  < 2e-16 ***
SoilTypeLow               0.178515   0.081896   2.180 0.029339 *  
SoilTypeMedium            0.215316   0.037505   5.741 1.02e-08 ***
LandTypeMedium Land       0.047731   0.027150   1.758 0.078826 .  
LandTypeUpland           -0.103589   0.055448  -1.868 0.061809 .  
as.factor(Year)2017-18   -0.235494   0.039795  -5.918 3.57e-09 ***
as.factor(Year)2018-19   -0.191806   0.037252  -5.149 2.76e-07 ***
as.factor(Year)2019-20   -0.130341   0.039874  -3.269 0.001090 ** 
as.factor(Year)2020-21   -0.221145   0.040089  -5.516 3.70e-08 ***
DistrictBegusarai        -0.013958   0.055308  -0.252 0.800768    
DistrictBuxar             0.069152   0.045503   1.520 0.128670    
DistrictDeoria           -0.475402   0.047403 -10.029  < 2e-16 ***
DistrictEast Champaran    0.081392   0.043326   1.879 0.060380 .  
DistrictKushinagar       -0.139145   0.050193  -2.772 0.005596 ** 
DistrictLakhisarai       -0.096428   0.053214  -1.812 0.070053 .  
DistrictMadhepura        -0.634577   0.046349 -13.691  < 2e-16 ***
DistrictMuzaffarpur      -0.021590   0.050784  -0.425 0.670758    
DistrictRohtas            0.610514   0.060737  10.052  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.6054 on 3621 degrees of freedom
Multiple R-squared:  0.6606,    Adjusted R-squared:  0.6582 
F-statistic: 271.1 on 26 and 3621 DF,  p-value: < 2.2e-16

# Breusch pagan test for heteroskedasticity
library(lmtest)
bptest(ols)

    studentized Breusch-Pagan test

data:  ols
BP = 303.78, df = 26, p-value < 2.2e-16

CSISA_KVKestim$ols_u=ols$residuals

CSISA_KVKestim$log_ols_u_sqrd=log((CSISA_KVKestim$ols_u)^2)

ols_log_variance = lm(log_ols_u_sqrd ~ SowingSchedule + VarietyClass_LDV + SeedRate + CropEstablishment + IrrigationNumber + SoilType + LandType + as.factor(Year) + District, data = CSISA_KVKestim)
summary(ols_log_variance)

Call:
lm(formula = log_ols_u_sqrd ~ SowingSchedule + VarietyClass_LDV + 
    SeedRate + CropEstablishment + IrrigationNumber + SoilType + 
    LandType + as.factor(Year) + District, data = CSISA_KVKestim)

Residuals:
     Min       1Q   Median       3Q      Max 
-14.1378  -1.0386   0.5213   1.4610   4.0038 

Coefficients:
                         Estimate Std. Error t value Pr(>|t|)    
(Intercept)              -2.28710    0.61550  -3.716 0.000206 ***
SowingScheduleT2          0.04563    0.14114   0.323 0.746492    
SowingScheduleT3         -0.02897    0.14602  -0.198 0.842725    
SowingScheduleT4          0.21919    0.15271   1.435 0.151286    
SowingScheduleT5          0.20446    0.16080   1.271 0.203642    
VarietyClass_LDV          0.06623    0.09415   0.703 0.481790    
SeedRate                  0.02936    0.01106   2.654 0.007993 ** 
CropEstablishmentCT-line -0.23767    0.64263  -0.370 0.711520    
CropEstablishmentZT      -0.93376    0.28507  -3.275 0.001065 ** 
IrrigationNumber         -0.03279    0.06807  -0.482 0.630079    
SoilTypeLow              -0.83666    0.30256  -2.765 0.005717 ** 
SoilTypeMedium           -0.30802    0.13856  -2.223 0.026279 *  
LandTypeMedium Land      -0.09958    0.10031  -0.993 0.320889    
LandTypeUpland           -0.07563    0.20485  -0.369 0.712020    
as.factor(Year)2017-18   -0.46565    0.14702  -3.167 0.001552 ** 
as.factor(Year)2018-19   -0.17665    0.13763  -1.284 0.199394    
as.factor(Year)2019-20   -0.37713    0.14731  -2.560 0.010506 *  
as.factor(Year)2020-21   -0.37383    0.14811  -2.524 0.011643 *  
DistrictBegusarai        -0.31749    0.20434  -1.554 0.120331    
DistrictBuxar             0.72204    0.16811   4.295 1.79e-05 ***
DistrictDeoria            0.13303    0.17513   0.760 0.447554    
DistrictEast Champaran   -0.23871    0.16007  -1.491 0.135972    
DistrictKushinagar        0.15601    0.18544   0.841 0.400237    
DistrictLakhisarai       -0.09794    0.19660  -0.498 0.618397    
DistrictMadhepura        -0.08121    0.17124  -0.474 0.635336    
DistrictMuzaffarpur       0.67695    0.18762   3.608 0.000313 ***
DistrictRohtas           -0.23035    0.22439  -1.027 0.304703    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.237 on 3621 degrees of freedom
Multiple R-squared:  0.03741,   Adjusted R-squared:  0.0305 
F-statistic: 5.412 on 26 and 3621 DF,  p-value: < 2.2e-16

# fitted values
CSISA_KVKestim$g_hat=fitted(ols_log_variance)
# Exponential of fitted values
CSISA_KVKestim$h_hat=exp(CSISA_KVKestim$g_hat)
# create weights
CSISA_KVKestim$w=1/CSISA_KVKestim$h_hat

# FGLS
fgls = lm(GrainYield ~ SowingSchedule +VarietyClass_LDV + SeedRate + CropEstablishment + IrrigationNumber + SoilType + LandType +as.factor(Year)+ District, data = CSISA_KVKestim, weight = w)
summary(fgls)

Call:
lm(formula = GrainYield ~ SowingSchedule + VarietyClass_LDV + 
    SeedRate + CropEstablishment + IrrigationNumber + SoilType + 
    LandType + as.factor(Year) + District, data = CSISA_KVKestim, 
    weights = w)

Weighted Residuals:
    Min      1Q  Median      3Q     Max 
-6.5729 -1.2736 -0.0168  1.2824  7.5161 

Coefficients:
                          Estimate Std. Error t value Pr(>|t|)    
(Intercept)               2.669334   0.181139  14.736  < 2e-16 ***
SowingScheduleT2         -0.294066   0.033419  -8.799  < 2e-16 ***
SowingScheduleT3         -0.627452   0.034741 -18.061  < 2e-16 ***
SowingScheduleT4         -1.056936   0.038154 -27.702  < 2e-16 ***
SowingScheduleT5         -1.376626   0.040011 -34.406  < 2e-16 ***
VarietyClass_LDV          0.435970   0.023947  18.206  < 2e-16 ***
SeedRate                 -0.006412   0.002751  -2.331  0.01981 *  
CropEstablishmentCT-line  0.924593   0.228001   4.055 5.11e-05 ***
CropEstablishmentZT       0.867418   0.118098   7.345 2.53e-13 ***
IrrigationNumber          0.455721   0.016756  27.197  < 2e-16 ***
SoilTypeLow               0.088969   0.067682   1.315  0.18876    
SoilTypeMedium            0.165849   0.037886   4.378 1.23e-05 ***
LandTypeMedium Land       0.034230   0.025991   1.317  0.18793    
LandTypeUpland           -0.131406   0.049275  -2.667  0.00769 ** 
as.factor(Year)2017-18   -0.116993   0.039639  -2.951  0.00318 ** 
as.factor(Year)2018-19   -0.116407   0.039095  -2.978  0.00292 ** 
as.factor(Year)2019-20    0.026882   0.040062   0.671  0.50226    
as.factor(Year)2020-21   -0.081225   0.040821  -1.990  0.04669 *  
DistrictBegusarai        -0.005822   0.047022  -0.124  0.90148    
DistrictBuxar             0.026466   0.049100   0.539  0.58991    
DistrictDeoria           -0.439702   0.043987  -9.996  < 2e-16 ***
DistrictEast Champaran    0.119128   0.038390   3.103  0.00193 ** 
DistrictKushinagar       -0.079281   0.046975  -1.688  0.09155 .  
DistrictLakhisarai       -0.146629   0.048949  -2.996  0.00276 ** 
DistrictMadhepura        -0.630738   0.041187 -15.314  < 2e-16 ***
DistrictMuzaffarpur      -0.012255   0.054433  -0.225  0.82188    
DistrictRohtas            0.611227   0.059065  10.348  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.89 on 3621 degrees of freedom
Multiple R-squared:  0.6972,    Adjusted R-squared:  0.695 
F-statistic: 320.6 on 26 and 3621 DF,  p-value: < 2.2e-16

library(stargazer)
stargazer(ols, fgls,ols_log_variance,
          column.labels=c("OLS","FGLS","ols_log_variance"),
          type="text",out ="JP_results.html",
          keep.stat=c("n","rsq"))

===========================================================
                                Dependent variable:        
                         ----------------------------------
                             GrainYield      log_ols_u_sqrd
                            OLS      FGLS         ols      
                            (1)       (2)         (3)      
-----------------------------------------------------------
SowingScheduleT2         -0.301*** -0.294***     0.046     
                          (0.038)   (0.033)     (0.141)    
                                                           
SowingScheduleT3         -0.640*** -0.627***     -0.029    
                          (0.040)   (0.035)     (0.146)    
                                                           
SowingScheduleT4         -1.039*** -1.057***     0.219     
                          (0.041)   (0.038)     (0.153)    
                                                           
SowingScheduleT5         -1.353*** -1.377***     0.204     
                          (0.044)   (0.040)     (0.161)    
                                                           
VarietyClass_LDV         0.441***  0.436***      0.066     
                          (0.025)   (0.024)     (0.094)    
                                                           
SeedRate                  -0.005*  -0.006**     0.029***   
                          (0.003)   (0.003)     (0.011)    
                                                           
CropEstablishmentCT-line 0.641***  0.925***      -0.238    
                          (0.174)   (0.228)     (0.643)    
                                                           
CropEstablishmentZT      0.663***  0.867***    -0.934***   
                          (0.077)   (0.118)     (0.285)    
                                                           
IrrigationNumber         0.493***  0.456***      -0.033    
                          (0.018)   (0.017)     (0.068)    
                                                           
SoilTypeLow               0.179**    0.089     -0.837***   
                          (0.082)   (0.068)     (0.303)    
                                                           
SoilTypeMedium           0.215***  0.166***     -0.308**   
                          (0.038)   (0.038)     (0.139)    
                                                           
LandTypeMedium Land       0.048*     0.034       -0.100    
                          (0.027)   (0.026)     (0.100)    
                                                           
LandTypeUpland            -0.104*  -0.131***     -0.076    
                          (0.055)   (0.049)     (0.205)    
                                                           
as.factor(Year)2017-18   -0.235*** -0.117***   -0.466***   
                          (0.040)   (0.040)     (0.147)    
                                                           
as.factor(Year)2018-19   -0.192*** -0.116***     -0.177    
                          (0.037)   (0.039)     (0.138)    
                                                           
as.factor(Year)2019-20   -0.130***   0.027      -0.377**   
                          (0.040)   (0.040)     (0.147)    
                                                           
as.factor(Year)2020-21   -0.221*** -0.081**     -0.374**   
                          (0.040)   (0.041)     (0.148)    
                                                           
DistrictBegusarai         -0.014    -0.006       -0.317    
                          (0.055)   (0.047)     (0.204)    
                                                           
DistrictBuxar              0.069     0.026      0.722***   
                          (0.046)   (0.049)     (0.168)    
                                                           
DistrictDeoria           -0.475*** -0.440***     0.133     
                          (0.047)   (0.044)     (0.175)    
                                                           
DistrictEast Champaran    0.081*   0.119***      -0.239    
                          (0.043)   (0.038)     (0.160)    
                                                           
DistrictKushinagar       -0.139***  -0.079*      0.156     
                          (0.050)   (0.047)     (0.185)    
                                                           
DistrictLakhisarai        -0.096*  -0.147***     -0.098    
                          (0.053)   (0.049)     (0.197)    
                                                           
DistrictMadhepura        -0.635*** -0.631***     -0.081    
                          (0.046)   (0.041)     (0.171)    
                                                           
DistrictMuzaffarpur       -0.022    -0.012      0.677***   
                          (0.051)   (0.054)     (0.188)    
                                                           
DistrictRohtas           0.611***  0.611***      -0.230    
                          (0.061)   (0.059)     (0.224)    
                                                           
Constant                 2.772***  2.669***    -2.287***   
                          (0.167)   (0.181)     (0.616)    
                                                           
-----------------------------------------------------------
Observations               3,648     3,648       3,648     
R2                         0.661     0.697       0.037     
===========================================================
Note:                           *p<0.1; **p<0.05; ***p<0.01

library(modelsummary)
list_models = list("Mean model:OLS"=ols,"Mean model: FGLS"=fgls, "Log Variance/Risk model"=ols_log_variance)

b <- list(geom_vline(xintercept = 0, color = 'orange'))
modelplot(list_models,background = b,coef_omit = "Interc")

modelsummary(list_models, stars=TRUE,output="tables/JP_results.docx")

2.1 Contribution of sowing dates and varieties to yield and risk

While the estimates shows the marginal effects of sowing dates or variety duration, one may be interest to know the contribution of these variables to the yield and yield risk. We use analysis of variance and shappley value regression. ### ANOVA

anova(ols)

Analysis of Variance Table

Response: GrainYield
                    Df  Sum Sq Mean Sq   F value    Pr(>F)    
SowingSchedule       4 1731.38  432.84 1180.8898 < 2.2e-16 ***
VarietyClass_LDV     1  160.32  160.32  437.3857 < 2.2e-16 ***
SeedRate             1   32.08   32.08   87.5170 < 2.2e-16 ***
CropEstablishment    2   12.40    6.20   16.9098 4.901e-08 ***
IrrigationNumber     1  360.23  360.23  982.7804 < 2.2e-16 ***
SoilType             2   11.95    5.97   16.2945 9.017e-08 ***
LandType             2   16.33    8.16   22.2729 2.432e-10 ***
as.factor(Year)      4    7.77    1.94    5.3013 0.0002958 ***
District             9  250.82   27.87   76.0322 < 2.2e-16 ***
Residuals         3621 1327.24    0.37                        
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

anova(fgls)

Analysis of Variance Table

Response: GrainYield
                    Df  Sum Sq Mean Sq  F value    Pr(>F)    
SowingSchedule       4 20279.6  5069.9 1419.573 < 2.2e-16 ***
VarietyClass_LDV     1  1645.2  1645.2  460.648 < 2.2e-16 ***
SeedRate             1   372.7   372.7  104.364 < 2.2e-16 ***
CropEstablishment    2    93.5    46.8   13.094 2.158e-06 ***
IrrigationNumber     1  3673.3  3673.3 1028.525 < 2.2e-16 ***
SoilType             2   288.1   144.1   40.339 < 2.2e-16 ***
LandType             2   240.6   120.3   33.690 3.186e-15 ***
as.factor(Year)      4    82.5    20.6    5.778  0.000124 ***
District             9  3096.3   344.0   96.330 < 2.2e-16 ***
Residuals         3621 12932.1     3.6                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

anova(ols_log_variance)

Analysis of Variance Table

Response: log_ols_u_sqrd
                    Df  Sum Sq Mean Sq F value    Pr(>F)    
SowingSchedule       4    50.7  12.684  2.5354 0.0382922 *  
VarietyClass_LDV     1     5.6   5.565  1.1123 0.2916550    
SeedRate             1     1.8   1.811  0.3620 0.5474131    
CropEstablishment    2   110.3  55.169 11.0272 1.681e-05 ***
IrrigationNumber     1     4.3   4.306  0.8607 0.3536113    
SoilType             2    38.1  19.049  3.8075 0.0222923 *  
LandType             2     3.9   1.942  0.3881 0.6783869    
as.factor(Year)      4   113.7  28.420  5.6806 0.0001482 ***
District             9   375.6  41.735  8.3420 2.093e-12 ***
Residuals         3621 18115.8   5.003                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

2.1.1 Shapley value regression

## Shapley value regression -----
library(ShapleyValue)

y <- CSISA_KVKestim$GrainYield
x=subset(CSISA_KVKestim, select=c("SowingSchedule","VarietyClass_LDV","SeedRate","CropEstablishment","IrrigationNumber","SoilType","LandType","Year","District"))

value <- shapleyvalue(y,x)

library(kableExtra)
value %>%
  kbl() %>%
  kable_classic(full_width = F, html_font = "Cambria")

	SowingSchedule	VarietyClass_LDV	SeedRate	CropEstablishment	IrrigationNumber	SoilType	LandType	Year	District
Shapley Value	0.2533	0.0944	0.0219	0.0054	0.1877	0.0042	0.0074	0.0063	0.0800
Standardized Shapley Value	0.3834	0.1429	0.0332	0.0082	0.2842	0.0064	0.0111	0.0095	0.1211

shapleyvaluet=as.data.frame(t(value))
shapleyvaluet=cbind(rownames(shapleyvaluet), data.frame(shapleyvaluet, row.names=NULL))
names(shapleyvaluet)[1]="vars"

library(ggplot2)
shapleyvalueplot=ggplot(shapleyvaluet,aes(x=reorder(vars,Standardized.Shapley.Value),y=Standardized.Shapley.Value))+
  geom_jitter(color="steelblue")+
  coord_flip()+
  labs(x="Variables",y="Standardized.Shapley.Value")
previous_theme <- theme_set(theme_bw())
shapleyvalueplot

ggsave("figures/Mean_shapleyvalueplot.png",dpi=300)


# Variance model
y <- CSISA_KVKestim$log_ols_u_sqrd
x=subset(CSISA_KVKestim, select=c("SowingSchedule","VarietyClass_LDV","SeedRate","CropEstablishment","IrrigationNumber","SoilType","LandType","Year","District"))

value <- shapleyvalue(y,x)

library(kableExtra)
value %>%
  kbl() %>%
  kable_classic(full_width = F, html_font = "Cambria")

	SowingSchedule	VarietyClass_LDV	SeedRate	CropEstablishment	IrrigationNumber	SoilType	LandType	Year	District
Shapley Value	0.0024	0.0001	0.0014	0.0044	0.0002	0.0029	0.0008	0.0060	0.0194
Standardized Shapley Value	0.0629	0.0028	0.0365	0.1168	0.0057	0.0765	0.0204	0.1611	0.5173

shapleyvaluet=as.data.frame(t(value))
shapleyvaluet=cbind(rownames(shapleyvaluet), data.frame(shapleyvaluet, row.names=NULL))
names(shapleyvaluet)[1]="vars"

library(ggplot2)
shapleyvalueplot=ggplot(shapleyvaluet,aes(x=reorder(vars,Standardized.Shapley.Value),y=Standardized.Shapley.Value))+
  geom_jitter(color="steelblue")+
  coord_flip()+
  labs(x="Variables",y="Standardized.Shapley.Value")
previous_theme <- theme_set(theme_bw())
shapleyvalueplot

ggsave("figures/Variance_shapleyvalueplot.png",dpi=300)

3 Moment approach

Antle (1983) extended the Just-Pope algorithm by including skewness as a measure of downside risk. The reason for the extension was the the Just-Pope approach considers variance as a measure of risk but variance doesnot distinguish unexpected bad events and good events. Skewness allows characterization of the unexpected downside effects.

ols_mean=lm(GrainYield ~ SowingSchedule+VarietyClass_LDV + SeedRate + CropEstablishment + IrrigationNumber + SoilType + LandType + as.factor(Year)+ District, data = CSISA_KVKestim)
summary(ols_mean)

Call:
lm(formula = GrainYield ~ SowingSchedule + VarietyClass_LDV + 
    SeedRate + CropEstablishment + IrrigationNumber + SoilType + 
    LandType + as.factor(Year) + District, data = CSISA_KVKestim)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.56769 -0.38434 -0.00054  0.40224  2.43981 

Coefficients:
                          Estimate Std. Error t value Pr(>|t|)    
(Intercept)               2.772121   0.166600  16.639  < 2e-16 ***
SowingScheduleT2         -0.300546   0.038204  -7.867 4.76e-15 ***
SowingScheduleT3         -0.640101   0.039523 -16.196  < 2e-16 ***
SowingScheduleT4         -1.038753   0.041335 -25.130  < 2e-16 ***
SowingScheduleT5         -1.353139   0.043525 -31.089  < 2e-16 ***
VarietyClass_LDV          0.441173   0.025483  17.313  < 2e-16 ***
SeedRate                 -0.005214   0.002994  -1.741 0.081713 .  
CropEstablishmentCT-line  0.641219   0.173942   3.686 0.000231 ***
CropEstablishmentZT       0.662518   0.077162   8.586  < 2e-16 ***
IrrigationNumber          0.492991   0.018425  26.756  < 2e-16 ***
SoilTypeLow               0.178515   0.081896   2.180 0.029339 *  
SoilTypeMedium            0.215316   0.037505   5.741 1.02e-08 ***
LandTypeMedium Land       0.047731   0.027150   1.758 0.078826 .  
LandTypeUpland           -0.103589   0.055448  -1.868 0.061809 .  
as.factor(Year)2017-18   -0.235494   0.039795  -5.918 3.57e-09 ***
as.factor(Year)2018-19   -0.191806   0.037252  -5.149 2.76e-07 ***
as.factor(Year)2019-20   -0.130341   0.039874  -3.269 0.001090 ** 
as.factor(Year)2020-21   -0.221145   0.040089  -5.516 3.70e-08 ***
DistrictBegusarai        -0.013958   0.055308  -0.252 0.800768    
DistrictBuxar             0.069152   0.045503   1.520 0.128670    
DistrictDeoria           -0.475402   0.047403 -10.029  < 2e-16 ***
DistrictEast Champaran    0.081392   0.043326   1.879 0.060380 .  
DistrictKushinagar       -0.139145   0.050193  -2.772 0.005596 ** 
DistrictLakhisarai       -0.096428   0.053214  -1.812 0.070053 .  
DistrictMadhepura        -0.634577   0.046349 -13.691  < 2e-16 ***
DistrictMuzaffarpur      -0.021590   0.050784  -0.425 0.670758    
DistrictRohtas            0.610514   0.060737  10.052  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.6054 on 3621 degrees of freedom
Multiple R-squared:  0.6606,    Adjusted R-squared:  0.6582 
F-statistic: 271.1 on 26 and 3621 DF,  p-value: < 2.2e-16

ols_variance=lm(ols_u^2~SowingSchedule+VarietyClass_LDV + SeedRate + CropEstablishment + IrrigationNumber + SoilType + LandType + as.factor(Year)+ District, data = CSISA_KVKestim)
summary(ols_variance)

Call:
lm(formula = ols_u^2 ~ SowingSchedule + VarietyClass_LDV + SeedRate + 
    CropEstablishment + IrrigationNumber + SoilType + LandType + 
    as.factor(Year) + District, data = CSISA_KVKestim)

Residuals:
    Min      1Q  Median      3Q     Max 
-0.7666 -0.2941 -0.1437  0.0914  5.8781 

Coefficients:
                          Estimate Std. Error t value Pr(>|t|)    
(Intercept)               0.026239   0.147010   0.178 0.858351    
SowingScheduleT2          0.059411   0.033711   1.762 0.078094 .  
SowingScheduleT3          0.012820   0.034876   0.368 0.713202    
SowingScheduleT4          0.128153   0.036475   3.513 0.000448 ***
SowingScheduleT5          0.189214   0.038407   4.927 8.75e-07 ***
VarietyClass_LDV         -0.047251   0.022486  -2.101 0.035681 *  
SeedRate                  0.007725   0.002642   2.924 0.003480 ** 
CropEstablishmentCT-line  0.002554   0.153490   0.017 0.986725    
CropEstablishmentZT      -0.125401   0.068089  -1.842 0.065599 .  
IrrigationNumber          0.047345   0.016259   2.912 0.003614 ** 
SoilTypeLow              -0.185661   0.072266  -2.569 0.010236 *  
SoilTypeMedium           -0.021887   0.033095  -0.661 0.508448    
LandTypeMedium Land      -0.059028   0.023958  -2.464 0.013794 *  
LandTypeUpland           -0.066748   0.048928  -1.364 0.172582    
as.factor(Year)2017-18   -0.152943   0.035116  -4.355 1.37e-05 ***
as.factor(Year)2018-19   -0.021560   0.032872  -0.656 0.511933    
as.factor(Year)2019-20   -0.050702   0.035186  -1.441 0.149674    
as.factor(Year)2020-21   -0.109369   0.035375  -3.092 0.002005 ** 
DistrictBegusarai        -0.045846   0.048805  -0.939 0.347606    
DistrictBuxar             0.346208   0.040153   8.622  < 2e-16 ***
DistrictDeoria            0.018077   0.041829   0.432 0.665652    
DistrictEast Champaran   -0.069652   0.038232  -1.822 0.068562 .  
DistrictKushinagar        0.020793   0.044291   0.469 0.638771    
DistrictLakhisarai        0.021967   0.046956   0.468 0.639944    
DistrictMadhepura         0.047363   0.040899   1.158 0.246924    
DistrictMuzaffarpur       0.215427   0.044813   4.807 1.59e-06 ***
DistrictRohtas            0.016187   0.053595   0.302 0.762647    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.5342 on 3621 degrees of freedom
Multiple R-squared:  0.08327,   Adjusted R-squared:  0.07669 
F-statistic: 12.65 on 26 and 3621 DF,  p-value: < 2.2e-16

ols_skewness=lm(ols_u^3~SowingSchedule+VarietyClass_LDV + SeedRate + CropEstablishment + IrrigationNumber + SoilType + LandType + as.factor(Year)+ District, data = CSISA_KVKestim)
summary(ols_skewness)

Call:
lm(formula = ols_u^3 ~ SowingSchedule + VarietyClass_LDV + SeedRate + 
    CropEstablishment + IrrigationNumber + SoilType + LandType + 
    as.factor(Year) + District, data = CSISA_KVKestim)

Residuals:
     Min       1Q   Median       3Q      Max 
-16.5786  -0.1467   0.0122   0.1718  14.4038 

Coefficients:
                           Estimate Std. Error t value Pr(>|t|)   
(Intercept)               1.331e-01  2.655e-01   0.501   0.6162   
SowingScheduleT2         -2.934e-02  6.088e-02  -0.482   0.6299   
SowingScheduleT3         -3.109e-02  6.299e-02  -0.494   0.6216   
SowingScheduleT4          1.034e-01  6.587e-02   1.570   0.1165   
SowingScheduleT5          1.231e-01  6.936e-02   1.774   0.0761 . 
VarietyClass_LDV         -2.227e-02  4.061e-02  -0.548   0.5834   
SeedRate                  1.539e-05  4.771e-03   0.003   0.9974   
CropEstablishmentCT-line -5.131e-01  2.772e-01  -1.851   0.0643 . 
CropEstablishmentZT      -3.956e-01  1.230e-01  -3.217   0.0013 **
IrrigationNumber          3.120e-02  2.936e-02   1.062   0.2881   
SoilTypeLow               1.854e-01  1.305e-01   1.421   0.1555   
SoilTypeMedium            1.303e-01  5.977e-02   2.179   0.0294 * 
LandTypeMedium Land       6.540e-02  4.327e-02   1.512   0.1307   
LandTypeUpland            1.084e-01  8.836e-02   1.227   0.2200   
as.factor(Year)2017-18    2.080e-02  6.342e-02   0.328   0.7429   
as.factor(Year)2018-19    9.287e-02  5.937e-02   1.564   0.1178   
as.factor(Year)2019-20   -1.369e-01  6.355e-02  -2.154   0.0313 * 
as.factor(Year)2020-21    8.989e-02  6.389e-02   1.407   0.1595   
DistrictBegusarai         1.338e-02  8.814e-02   0.152   0.8793   
DistrictBuxar            -1.238e-01  7.252e-02  -1.708   0.0878 . 
DistrictDeoria           -8.820e-02  7.554e-02  -1.168   0.2431   
DistrictEast Champaran    2.585e-03  6.905e-02   0.037   0.9701   
DistrictKushinagar       -9.218e-02  7.999e-02  -1.152   0.2492   
DistrictLakhisarai        4.759e-02  8.480e-02   0.561   0.5747   
DistrictMadhepura        -3.641e-03  7.386e-02  -0.049   0.9607   
DistrictMuzaffarpur      -1.022e-01  8.093e-02  -1.262   0.2069   
DistrictRohtas            1.814e-03  9.679e-02   0.019   0.9851   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.9648 on 3621 degrees of freedom
Multiple R-squared:  0.01842,   Adjusted R-squared:  0.01137 
F-statistic: 2.613 on 26 and 3621 DF,  p-value: 1.508e-05

ols_kurtosis=lm(ols_u^4~SowingSchedule+VarietyClass_LDV + SeedRate + CropEstablishment + IrrigationNumber + SoilType + LandType + as.factor(Year)+ District, data = CSISA_KVKestim)
summary(ols_kurtosis)

Call:
lm(formula = ols_u^4 ~ SowingSchedule + VarietyClass_LDV + SeedRate + 
    CropEstablishment + IrrigationNumber + SoilType + LandType + 
    as.factor(Year) + District, data = CSISA_KVKestim)

Residuals:
   Min     1Q Median     3Q    Max 
-1.890 -0.519 -0.247  0.044 41.895 

Coefficients:
                         Estimate Std. Error t value Pr(>|t|)    
(Intercept)              -0.78854    0.43850  -1.798 0.072219 .  
SowingScheduleT2          0.21106    0.10055   2.099 0.035886 *  
SowingScheduleT3          0.06119    0.10403   0.588 0.556407    
SowingScheduleT4          0.28845    0.10880   2.651 0.008053 ** 
SowingScheduleT5          0.48739    0.11456   4.254 2.15e-05 ***
VarietyClass_LDV         -0.22727    0.06707  -3.388 0.000710 ***
SeedRate                  0.01823    0.00788   2.313 0.020793 *  
CropEstablishmentCT-line  0.08120    0.45783   0.177 0.859240    
CropEstablishmentZT      -0.05976    0.20309  -0.294 0.768599    
IrrigationNumber          0.15189    0.04850   3.132 0.001750 ** 
SoilTypeLow              -0.29516    0.21555  -1.369 0.170997    
SoilTypeMedium            0.05895    0.09871   0.597 0.550410    
LandTypeMedium Land      -0.16283    0.07146  -2.279 0.022754 *  
LandTypeUpland           -0.20286    0.14594  -1.390 0.164619    
as.factor(Year)2017-18   -0.38553    0.10474  -3.681 0.000236 ***
as.factor(Year)2018-19   -0.06095    0.09805  -0.622 0.534252    
as.factor(Year)2019-20    0.02559    0.10495   0.244 0.807372    
as.factor(Year)2020-21   -0.26246    0.10552  -2.487 0.012913 *  
DistrictBegusarai        -0.02043    0.14557  -0.140 0.888421    
DistrictBuxar             0.98589    0.11977   8.232 2.55e-16 ***
DistrictDeoria            0.01428    0.12477   0.114 0.908867    
DistrictEast Champaran   -0.06788    0.11404  -0.595 0.551718    
DistrictKushinagar        0.07752    0.13211   0.587 0.557389    
DistrictLakhisarai        0.07637    0.14006   0.545 0.585616    
DistrictMadhepura         0.15332    0.12199   1.257 0.208895    
DistrictMuzaffarpur       0.49656    0.13367   3.715 0.000206 ***
DistrictRohtas            0.16020    0.15986   1.002 0.316340    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.594 on 3621 degrees of freedom
Multiple R-squared:  0.06319,   Adjusted R-squared:  0.05647 
F-statistic: 9.395 on 26 and 3621 DF,  p-value: < 2.2e-16

library(stargazer)
stargazer(ols_mean, ols_variance,ols_skewness,ols_kurtosis,
          column.labels=c("Mean","Variance","Skewness","Kurtosis"),
          type="text", out ="Moments_results.html",
          keep.stat=c("n","rsq"))

=================================================================
                                   Dependent variable:           
                         ----------------------------------------
                         GrainYield  ols_u2    ols_u3    ols_u4  
                            Mean    Variance  Skewness  Kurtosis 
                            (1)        (2)       (3)       (4)   
-----------------------------------------------------------------
SowingScheduleT2         -0.301***   0.059*    -0.029    0.211** 
                          (0.038)    (0.034)   (0.061)   (0.101) 
                                                                 
SowingScheduleT3         -0.640***    0.013    -0.031     0.061  
                          (0.040)    (0.035)   (0.063)   (0.104) 
                                                                 
SowingScheduleT4         -1.039***  0.128***    0.103   0.288*** 
                          (0.041)    (0.036)   (0.066)   (0.109) 
                                                                 
SowingScheduleT5         -1.353***  0.189***   0.123*   0.487*** 
                          (0.044)    (0.038)   (0.069)   (0.115) 
                                                                 
VarietyClass_LDV          0.441***  -0.047**   -0.022   -0.227***
                          (0.025)    (0.022)   (0.041)   (0.067) 
                                                                 
SeedRate                  -0.005*   0.008***   0.00002   0.018** 
                          (0.003)    (0.003)   (0.005)   (0.008) 
                                                                 
CropEstablishmentCT-line  0.641***    0.003    -0.513*    0.081  
                          (0.174)    (0.153)   (0.277)   (0.458) 
                                                                 
CropEstablishmentZT       0.663***   -0.125*  -0.396***  -0.060  
                          (0.077)    (0.068)   (0.123)   (0.203) 
                                                                 
IrrigationNumber          0.493***  0.047***    0.031   0.152*** 
                          (0.018)    (0.016)   (0.029)   (0.048) 
                                                                 
SoilTypeLow               0.179**   -0.186**    0.185    -0.295  
                          (0.082)    (0.072)   (0.131)   (0.216) 
                                                                 
SoilTypeMedium            0.215***   -0.022    0.130**    0.059  
                          (0.038)    (0.033)   (0.060)   (0.099) 
                                                                 
LandTypeMedium Land        0.048*   -0.059**    0.065   -0.163** 
                          (0.027)    (0.024)   (0.043)   (0.071) 
                                                                 
LandTypeUpland            -0.104*    -0.067     0.108    -0.203  
                          (0.055)    (0.049)   (0.088)   (0.146) 
                                                                 
as.factor(Year)2017-18   -0.235***  -0.153***   0.021   -0.386***
                          (0.040)    (0.035)   (0.063)   (0.105) 
                                                                 
as.factor(Year)2018-19   -0.192***   -0.022     0.093    -0.061  
                          (0.037)    (0.033)   (0.059)   (0.098) 
                                                                 
as.factor(Year)2019-20   -0.130***   -0.051   -0.137**    0.026  
                          (0.040)    (0.035)   (0.064)   (0.105) 
                                                                 
as.factor(Year)2020-21   -0.221***  -0.109***   0.090   -0.262** 
                          (0.040)    (0.035)   (0.064)   (0.106) 
                                                                 
DistrictBegusarai          -0.014    -0.046     0.013    -0.020  
                          (0.055)    (0.049)   (0.088)   (0.146) 
                                                                 
DistrictBuxar              0.069    0.346***   -0.124*  0.986*** 
                          (0.046)    (0.040)   (0.073)   (0.120) 
                                                                 
DistrictDeoria           -0.475***    0.018    -0.088     0.014  
                          (0.047)    (0.042)   (0.076)   (0.125) 
                                                                 
DistrictEast Champaran     0.081*    -0.070*    0.003    -0.068  
                          (0.043)    (0.038)   (0.069)   (0.114) 
                                                                 
DistrictKushinagar       -0.139***    0.021    -0.092     0.078  
                          (0.050)    (0.044)   (0.080)   (0.132) 
                                                                 
DistrictLakhisarai        -0.096*     0.022     0.048     0.076  
                          (0.053)    (0.047)   (0.085)   (0.140) 
                                                                 
DistrictMadhepura        -0.635***    0.047    -0.004     0.153  
                          (0.046)    (0.041)   (0.074)   (0.122) 
                                                                 
DistrictMuzaffarpur        -0.022   0.215***   -0.102   0.497*** 
                          (0.051)    (0.045)   (0.081)   (0.134) 
                                                                 
DistrictRohtas            0.611***    0.016     0.002     0.160  
                          (0.061)    (0.054)   (0.097)   (0.160) 
                                                                 
Constant                  2.772***    0.026     0.133    -0.789* 
                          (0.167)    (0.147)   (0.266)   (0.438) 
                                                                 
-----------------------------------------------------------------
Observations               3,648      3,648     3,648     3,648  
R2                         0.661      0.083     0.018     0.063  
=================================================================
Note:                                 *p<0.1; **p<0.05; ***p<0.01

library(modelsummary)
list_models_moments = list("Mean model"=ols_mean, "Variance model"=ols_variance,"Skewness model"=ols_skewness,"Kurtosis model"=ols_kurtosis)

modelsummary(list_models_moments,stars=TRUE)

	Mean model	Variance model	Skewness model	Kurtosis model
(Intercept)	2.772***	0.026	0.133	-0.789+
	(0.167)	(0.147)	(0.266)	(0.438)
SowingScheduleT2	-0.301***	0.059+	-0.029	0.211*
	(0.038)	(0.034)	(0.061)	(0.101)
SowingScheduleT3	-0.640***	0.013	-0.031	0.061
	(0.040)	(0.035)	(0.063)	(0.104)
SowingScheduleT4	-1.039***	0.128***	0.103	0.288**
	(0.041)	(0.036)	(0.066)	(0.109)
SowingScheduleT5	-1.353***	0.189***	0.123+	0.487***
	(0.044)	(0.038)	(0.069)	(0.115)
VarietyClass_LDV	0.441***	-0.047*	-0.022	-0.227***
	(0.025)	(0.022)	(0.041)	(0.067)
SeedRate	-0.005+	0.008**	0.000	0.018*
	(0.003)	(0.003)	(0.005)	(0.008)
CropEstablishmentCT-line	0.641***	0.003	-0.513+	0.081
	(0.174)	(0.153)	(0.277)	(0.458)
CropEstablishmentZT	0.663***	-0.125+	-0.396**	-0.060
	(0.077)	(0.068)	(0.123)	(0.203)
IrrigationNumber	0.493***	0.047**	0.031	0.152**
	(0.018)	(0.016)	(0.029)	(0.048)
SoilTypeLow	0.179*	-0.186*	0.185	-0.295
	(0.082)	(0.072)	(0.131)	(0.216)
SoilTypeMedium	0.215***	-0.022	0.130*	0.059
	(0.038)	(0.033)	(0.060)	(0.099)
LandTypeMedium Land	0.048+	-0.059*	0.065	-0.163*
	(0.027)	(0.024)	(0.043)	(0.071)
LandTypeUpland	-0.104+	-0.067	0.108	-0.203
	(0.055)	(0.049)	(0.088)	(0.146)
as.factor(Year)2017-18	-0.235***	-0.153***	0.021	-0.386***
	(0.040)	(0.035)	(0.063)	(0.105)
as.factor(Year)2018-19	-0.192***	-0.022	0.093	-0.061
	(0.037)	(0.033)	(0.059)	(0.098)
as.factor(Year)2019-20	-0.130**	-0.051	-0.137*	0.026
	(0.040)	(0.035)	(0.064)	(0.105)
as.factor(Year)2020-21	-0.221***	-0.109**	0.090	-0.262*
	(0.040)	(0.035)	(0.064)	(0.106)
DistrictBegusarai	-0.014	-0.046	0.013	-0.020
	(0.055)	(0.049)	(0.088)	(0.146)
DistrictBuxar	0.069	0.346***	-0.124+	0.986***
	(0.046)	(0.040)	(0.073)	(0.120)
DistrictDeoria	-0.475***	0.018	-0.088	0.014
	(0.047)	(0.042)	(0.076)	(0.125)
DistrictEast Champaran	0.081+	-0.070+	0.003	-0.068
	(0.043)	(0.038)	(0.069)	(0.114)
DistrictKushinagar	-0.139**	0.021	-0.092	0.078
	(0.050)	(0.044)	(0.080)	(0.132)
DistrictLakhisarai	-0.096+	0.022	0.048	0.076
	(0.053)	(0.047)	(0.085)	(0.140)
DistrictMadhepura	-0.635***	0.047	-0.004	0.153
	(0.046)	(0.041)	(0.074)	(0.122)
DistrictMuzaffarpur	-0.022	0.215***	-0.102	0.497***
	(0.051)	(0.045)	(0.081)	(0.134)
DistrictRohtas	0.611***	0.016	0.002	0.160
	(0.061)	(0.054)	(0.097)	(0.160)
Num.Obs.	3648	3648	3648	3648
R2	0.661	0.083	0.018	0.063
R2 Adj.	0.658	0.077	0.011	0.056
AIC	6720.2	5807.5	10120.3	13781.0
BIC	6893.8	5981.2	10294.0	13954.6
Log.Lik.	-3332.089	-2875.757	-5032.155	-6862.493
F	271.065	12.651	2.613	9.395
RMSE	0.60	0.53	0.96	1.59
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001

b <- list(geom_vline(xintercept = 0, color = 'orange'))

modelplot(list_models_moments,stars=TRUE, coef_omit = "Interc", background=b)

modelsummary(list_models_moments, stars=TRUE,output="tables/Moments_results.docx")

3.1 Contribution to moments

3.1.1 ANOVA

anova(ols_mean)

Analysis of Variance Table

Response: GrainYield
                    Df  Sum Sq Mean Sq   F value    Pr(>F)    
SowingSchedule       4 1731.38  432.84 1180.8898 < 2.2e-16 ***
VarietyClass_LDV     1  160.32  160.32  437.3857 < 2.2e-16 ***
SeedRate             1   32.08   32.08   87.5170 < 2.2e-16 ***
CropEstablishment    2   12.40    6.20   16.9098 4.901e-08 ***
IrrigationNumber     1  360.23  360.23  982.7804 < 2.2e-16 ***
SoilType             2   11.95    5.97   16.2945 9.017e-08 ***
LandType             2   16.33    8.16   22.2729 2.432e-10 ***
as.factor(Year)      4    7.77    1.94    5.3013 0.0002958 ***
District             9  250.82   27.87   76.0322 < 2.2e-16 ***
Residuals         3621 1327.24    0.37                        
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

anova(ols_variance)

Analysis of Variance Table

Response: ols_u^2
                    Df  Sum Sq Mean Sq F value    Pr(>F)    
SowingSchedule       4   17.04  4.2591 14.9227 4.230e-12 ***
VarietyClass_LDV     1    0.10  0.0988  0.3461  0.556380    
SeedRate             1    0.00  0.0043  0.0152  0.901902    
CropEstablishment    2    2.89  1.4428  5.0551  0.006422 ** 
IrrigationNumber     1    7.60  7.5970 26.6178 2.613e-07 ***
SoilType             2    1.47  0.7354  2.5767  0.076165 .  
LandType             2    1.72  0.8581  3.0064  0.049592 *  
as.factor(Year)      4   10.52  2.6297  9.2137 2.098e-07 ***
District             9   52.55  5.8389 20.4578 < 2.2e-16 ***
Residuals         3621 1033.47  0.2854                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

anova(ols_skewness)

Analysis of Variance Table

Response: ols_u^3
                    Df Sum Sq Mean Sq F value    Pr(>F)    
SowingSchedule       4   11.1  2.7831  2.9897   0.01779 *  
VarietyClass_LDV     1    0.0  0.0047  0.0050   0.94340    
SeedRate             1    0.8  0.8466  0.9095   0.34032    
CropEstablishment    2    6.1  3.0521  3.2786   0.03779 *  
IrrigationNumber     1    0.3  0.2568  0.2759   0.59943    
SoilType             2    7.4  3.6933  3.9674   0.01900 *  
LandType             2    3.0  1.4780  1.5876   0.20455    
as.factor(Year)      4   24.9  6.2304  6.6928 2.307e-05 ***
District             9    9.6  1.0705  1.1500   0.32328    
Residuals         3621 3370.8  0.9309                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

anova(ols_kurtosis)

Analysis of Variance Table

Response: ols_u^4
                    Df Sum Sq Mean Sq F value    Pr(>F)    
SowingSchedule       4  104.7  26.163 10.3034 2.698e-08 ***
VarietyClass_LDV     1   10.0   9.966  3.9249   0.04765 *  
SeedRate             1    0.0   0.013  0.0050   0.94354    
CropEstablishment    2    3.4   1.696  0.6681   0.51276    
IrrigationNumber     1   57.4  57.358 22.5884 2.085e-06 ***
SoilType             2    6.8   3.393  1.3363   0.26296    
LandType             2   16.7   8.356  3.2908   0.03733 *  
as.factor(Year)      4   72.4  18.101  7.1285 1.030e-05 ***
District             9  349.0  38.774 15.2696 < 2.2e-16 ***
Residuals         3621 9194.8   2.539                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

3.1.2 Shapley value regression

## Shapley value regression -----
library(ShapleyValue)

y <- CSISA_KVKestim$GrainYield
x=subset(CSISA_KVKestim, select=c("SowingSchedule","VarietyClass_LDV","SeedRate","CropEstablishment","IrrigationNumber","SoilType","LandType","Year","District"))

value <- shapleyvalue(y,x)

library(kableExtra)
value %>%
  kbl() %>%
  kable_classic(full_width = F, html_font = "Cambria")

	SowingSchedule	VarietyClass_LDV	SeedRate	CropEstablishment	IrrigationNumber	SoilType	LandType	Year	District
Shapley Value	0.2533	0.0944	0.0219	0.0054	0.1877	0.0042	0.0074	0.0063	0.0800
Standardized Shapley Value	0.3834	0.1429	0.0332	0.0082	0.2842	0.0064	0.0111	0.0095	0.1211

shapleyvaluet=as.data.frame(t(value))
shapleyvaluet=cbind(rownames(shapleyvaluet), data.frame(shapleyvaluet, row.names=NULL))
names(shapleyvaluet)[1]="vars"

library(ggplot2)
shapleyvalueplot=ggplot(shapleyvaluet,aes(x=reorder(vars,Standardized.Shapley.Value),y=Standardized.Shapley.Value))+
  geom_jitter(color="steelblue")+
  coord_flip()+
  labs(x="Variables",y="Standardized.Shapley.Value")
previous_theme <- theme_set(theme_bw())
shapleyvalueplot

ggsave("figures/Mean_shapleyvalueplot.png",dpi=300)


# Variance model
y <- CSISA_KVKestim$ols_u^2
x=subset(CSISA_KVKestim, select=c("SowingSchedule","VarietyClass_LDV","SeedRate","CropEstablishment","IrrigationNumber","SoilType","LandType","Year","District"))

value <- shapleyvalue(y,x)

library(kableExtra)
value %>%
  kbl() %>%
  kable_classic(full_width = F, html_font = "Cambria")

	SowingSchedule	VarietyClass_LDV	SeedRate	CropEstablishment	IrrigationNumber	SoilType	LandType	Year	District
Shapley Value	0.0140	0.0016	0.0017	0.0017	0.0019	0.0021	0.0030	0.0095	0.0478
Standardized Shapley Value	0.1684	0.0188	0.0205	0.0198	0.0230	0.0256	0.0357	0.1141	0.5740

shapleyvaluet=as.data.frame(t(value))
shapleyvaluet=cbind(rownames(shapleyvaluet), data.frame(shapleyvaluet, row.names=NULL))
names(shapleyvaluet)[1]="vars"

library(ggplot2)
shapleyvalueplot=ggplot(shapleyvaluet,aes(x=reorder(vars,Standardized.Shapley.Value),y=Standardized.Shapley.Value))+
  geom_jitter(color="steelblue")+
  coord_flip()+
  labs(x="Variables",y="Standardized.Shapley.Value")
previous_theme <- theme_set(theme_bw())
shapleyvalueplot

ggsave("figures/Variance_shapleyvalueplot.png",dpi=300)

# Skewness model
y <- CSISA_KVKestim$ols_u^3
x=subset(CSISA_KVKestim, select=c("SowingSchedule","VarietyClass_LDV","SeedRate","CropEstablishment","IrrigationNumber","SoilType","LandType","Year","District"))

value <- shapleyvalue(y,x)

library(kableExtra)
value %>%
  kbl() %>%
  kable_classic(full_width = F, html_font = "Cambria")

	SowingSchedule	VarietyClass_LDV	SeedRate	CropEstablishment	IrrigationNumber	SoilType	LandType	Year	District
Shapley Value	0.0034	0.0002	0.0003	0.0023	0.0002	0.0014	0.001	0.0072	0.0025
Standardized Shapley Value	0.1863	0.0089	0.0148	0.1269	0.0088	0.0737	0.055	0.3920	0.1336

shapleyvaluet=as.data.frame(t(value))
shapleyvaluet=cbind(rownames(shapleyvaluet), data.frame(shapleyvaluet, row.names=NULL))
names(shapleyvaluet)[1]="vars"

library(ggplot2)
shapleyvalueplot=ggplot(shapleyvaluet,aes(x=reorder(vars,Standardized.Shapley.Value),y=Standardized.Shapley.Value))+
  geom_jitter(color="steelblue")+
  coord_flip()+
  labs(x="Variables",y="Standardized.Shapley.Value")
previous_theme <- theme_set(theme_bw())
shapleyvalueplot

ggsave("figures/Skewness_shapleyvalueplot.png",dpi=300)

# Kurtosis model
y <- CSISA_KVKestim$ols_u^4
x=subset(CSISA_KVKestim, select=c("SowingSchedule","VarietyClass_LDV","SeedRate","CropEstablishment","IrrigationNumber","SoilType","LandType","Year","District"))

value <- shapleyvalue(y,x)

library(kableExtra)
value %>%
  kbl() %>%
  kable_classic(full_width = F, html_font = "Cambria")

	SowingSchedule	VarietyClass_LDV	SeedRate	CropEstablishment	IrrigationNumber	SoilType	LandType	Year	District
Shapley Value	0.0095	0.0032	0.0011	0.0002	0.0020	0.0009	0.0023	0.0079	0.0362
Standardized Shapley Value	0.1503	0.0503	0.0168	0.0032	0.0324	0.0136	0.0359	0.1242	0.5732

shapleyvaluet=as.data.frame(t(value))
shapleyvaluet=cbind(rownames(shapleyvaluet), data.frame(shapleyvaluet, row.names=NULL))
names(shapleyvaluet)[1]="vars"

library(ggplot2)
shapleyvalueplot=ggplot(shapleyvaluet,aes(x=reorder(vars,Standardized.Shapley.Value),y=Standardized.Shapley.Value))+
  geom_jitter(color="steelblue")+
  coord_flip()+
  labs(x="Variables",y="Standardized.Shapley.Value")
previous_theme <- theme_set(theme_bw())
shapleyvalueplot

ggsave("figures/Kurtosis_shapleyvalueplot.png",dpi=300)

4 Addressing endogeneity using Lewbel (2012) method

When using observational data, there may be endogeneity associated with the choice of the technology.

library(REndo)

hetEr <- hetErrorsIV(GrainYield~SowingSchedule_T1+SowingSchedule_T2+SowingSchedule_T3+SowingSchedule_T4+IrrigationNumber + SeedRate+VarietyClass_LDV |SowingSchedule_T1| IIV(SowingSchedule_T2,SowingSchedule_T3,SowingSchedule_T4,IrrigationNumber,SeedRate,VarietyClass_LDV ), data = CSISA_KVKestim)

summary(hetEr)

Call:
hetErrorsIV(formula = GrainYield ~ SowingSchedule_T1 + SowingSchedule_T2 + 
    SowingSchedule_T3 + SowingSchedule_T4 + IrrigationNumber + 
    SeedRate + VarietyClass_LDV | SowingSchedule_T1 | IIV(SowingSchedule_T2, 
    SowingSchedule_T3, SowingSchedule_T4, IrrigationNumber, SeedRate, 
    VarietyClass_LDV), data = CSISA_KVKestim)

Residuals:
      Min        1Q    Median        3Q       Max 
-2.316783 -0.439010  0.003074  0.459840  2.635782 

Coefficients:
                   Estimate Std. Error t value Pr(>|t|)    
(Intercept)        2.203195   0.153986  14.308  < 2e-16 ***
SowingSchedule_T1  1.328627   0.048057  27.647  < 2e-16 ***
SowingSchedule_T2  1.056724   0.038467  27.471  < 2e-16 ***
SowingSchedule_T3  0.728336   0.035591  20.464  < 2e-16 ***
SowingSchedule_T4  0.309741   0.034204   9.056  < 2e-16 ***
IrrigationNumber   0.541894   0.019000  28.521  < 2e-16 ***
SeedRate          -0.011693   0.002974  -3.931  8.6e-05 ***
VarietyClass_LDV   0.395841   0.027046  14.636  < 2e-16 ***

Diagnostic tests:
                  df1  df2 statistic  p-value    
Weak instruments    6 3635  6803.309  < 2e-16 ***
Wu-Hausman          1 3639     2.566    0.109    
Sargan              5   NA    54.066 2.03e-10 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.6669 on 3640 degrees of freedom
Multiple R-Squared: 0.586,  Adjusted R-squared: 0.5852 
Wald test: 730.1 on 7 and 3640 DF,  p-value: < 2.2e-16 

library(modelsummary)

modelplot(hetEr)

5 Spatial JP and Moments models

Risk exposure is usually spatially dependent. Farms close together in space are more likely to face simular exposures than farmers away. In addition, one may want to predict out of sample beyond locations where the agronomic trials were conducted. We show how to use spatially varying coefficient models to develop a surface of the effectiveness of each of the sowing date strategies and variety classes across Bihar. This approach allows us to recommend, with associated measures of uncertainty sowing dates and variety classes that maximize yieCSISA_KVKestim, minimize variance, skewness and kurtosis.

# Spatially varying coefficient model -----------------------------------------
# library(spBayes)
# coords=dplyr::select(CSISA_KVKestim,Longitude,Latitude)
# coords=as.matrix(coords)

# library(geoR)
# coords=jitterDupCoords(coords,min=2,max=10)

# nrow(coords)
# nrow(tau.hat_weeding)

# n.samples <- 1000 

# t1 <- Sys.time()

# r <-6
# n.ltr <- r*(r+1)/2

# priors <- list("phi.Unif"=list(rep(1,r), rep(10,r)), "K.IW"=list(r, diag(rep(1,r))), "tau.sq.IG"=c(2, 1))

# starting <- list("phi"=rep(3/0.5,r), "A"=rep(1,n.ltr), "tau.sq"=1) 

# tuning <- list("phi"=rep(0.1,r), "A"=rep(0.01, n.ltr), "tau.sq"=0.01)


# cf.weeding.spVC2 <- spBayes::spSVC(GrainYield~SowingSchedule_T1+SowingSchedule_T2+SowingSchedule_T3+SowingSchedule_T4+VarietyClass_LDV+IrrigationNumber + SeedRate, data=CSISA_KVKestim,coords=coords,
#                                   starting= starting,svc.cols=c(1,2,3,4,5,6),
#                                   tuning=tuning,
#                                   priors=priors,
#                                   cov.model="exponential",n.samples=n.samples,
#                                   n.omp.threads=15)


# t2 <- Sys.time()
# t2 - t1

# library(coda)
# round(summary(mcmc(cf.weeding.spVC2 $p.theta.samples))$quantiles,3)

# burn.in <- floor(0.75*n.samples) 

# cf.mean.sp.r <- spRecover(cf.weeding.spVC2 , start=burn.in,n.omp.threads=15)

# round(summary(cf.mean.sp.r$p.beta.recover.samples)$quantiles[,c(3,1,5)],2)

# round(summary(cf.mean.sp.r$p.theta.recover.samples)$quantiles[,c(3,1,5)],2)

# tilde.beta.0=apply(cf.mean.sp.r$p.tilde.beta.recover.samples[["tilde.beta.(Intercept)"]],1,median)

# tilde.beta.t1=apply(cf.mean.sp.r$p.tilde.beta.recover.samples[["tilde.beta.t1"]],1,median)



# library(terra)
# library(stars)
# library(raster)

# CSISA_KVKestim_sp= SpatialPointsDataFrame(cbind(CSISA_KVKestim$Longitude,CSISA_KVKestim$Latitude),data=CSISA_KVKestim,proj4string=CRS("+proj=longlat +datum=WGS84"))

# plot(CSISA_KVKestim_sp)

# library(geodata)
# India=gadm(country="IND", level=1, path="temp")
# plot(India)
# India_Bihar=subset(India,India$NAME_1=="Bihar")
# plot(India_Bihar)

# library(sf)
# India_Bihar=st_as_sf(India_Bihar)
# India_Bihar=as_Spatial(India_Bihar)
# # Inverse distance approach -----
# library(gstat) # Use gstat's idw routine
# library(sp)    # Used for the spsample function
# library(tmap)

# wgs84.prj=CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")
# India_Bihar_wgs84.prj= spTransform(India_Bihar,wgs84.prj)

# CSISA_KVKestim_sp@bbox <- India_Bihar_wgs84.prj@bbox 

# grd <- as.data.frame(spsample(CSISA_KVKestim_sp, "regular", n=10000))

# names(grd)       <- c("X", "Y")
# coordinates(grd) <- c("X", "Y")
# gridded(grd)     <- TRUE  # Create SpatialPixel object
# fullgrid(grd)    <- TRUE  # Create SpatialGrid object
# plot(grd)
# proj4string(CSISA_KVKestim_sp) <- proj4string(CSISA_KVKestim_sp) # Temp fix until new proj env is adopted
# proj4string(grd) <- proj4string(CSISA_KVKestim_sp)

# library(spBayes)

# India_Bihar_wgs84.prjpoly <- India_Bihar_wgs84.prj@polygons[[1]]@Polygons[[1]]@coords 
# India_Bihar_wgs84.prjpoly <- as.matrix(India_Bihar_wgs84.prjpoly)

# pred.coords <- SpatialPoints(grd)@coords 

# pred.covars <- as.matrix(rep(1, nrow(pred.coords)))

# cf.mean.sp.pred <- spPredict(cf.mean.sp.r,pred.coords=pred.coords,
#                                     pred.covars=pred.covars,n.omp.threads=15, thin=25, joint=TRUE)

# cf.mean.sp.pred.pred.mu = apply(cf.mean.sp.pred$p.y.predictive.samples,1,mean)
# cf.mean.sp.pred.sd = apply(cf.mean.sp.pred$p.y.predictive.samples,1,sd)

# x.res=100
# library(MBA)
# surf <- mba.surf(cbind(coords, tau.hat.cf_irrign), no.X=x.res, no.Y=x.res,extend=TRUE, sp=TRUE)$xyz.est
# surf <- as.image.SpatialGridDataFrame(surf) 
# z.lim <- range(surf[["z"]], na.rm=TRUE) 
# pred.grid <- as.data.frame(list(pred.coords,pred.mu=cf.mean.sp.pred.pred.mu,pred.sd=cf.mean.sp.pred.sd))

# coordinates(pred.grid) = c("X", "Y") 
# gridded(pred.grid) <- TRUE 
# pred.mu.image <- as.image.SpatialGridDataFrame(pred.grid["pred.mu"])
# pred.sd.image <- as.image.SpatialGridDataFrame(pred.grid["pred.sd"])

# library(fieCSISA_KVKestim)
# image.plot(surf, axes=TRUE, zlim=z.lim, col=tim.colors(25),xaxs = "r", yaxs = "r",main="Yield gain to two means")
# plot(India_State_Boundary_Bihar_wgs84.prj, add=TRUE) 



# png("figures/BayesianKriggedMeanValue_irrign2_mu.png")
# #par(mfrow=c(1,2)) 
# image.plot(pred.mu.image,xaxs = "r", yaxs = "r",main="Mean predicted yield gain to two means")
# plot(India_State_Boundary_Bihar_wgs84.prj, add=TRUE)
# dev.off()

# png("figures/BayesianKriggedMeanValue_irrign2_sd.png")
# image.plot(pred.sd.image,xaxs = "r", yaxs = "r",main="Sd of predicted yield gain to two means")
# plot(India_State_Boundary_Bihar_wgs84.prj, add=TRUE)
# dev.off()

# #pred.mu.image=crop(pred.mu.image,India_State_Boundary_Bihar_wgs84.prj)

# writeGDAL(pred.grid["pred.mu"], "figures/mean2.tau.pred.mu.image.tif") 
# writeGDAL(pred.grid["pred.sd"], "figures/mean2.tau.pred.sd.image.tif")


# library(rasterVis)
# pred.mu_irrign2=pred.grid["pred.mu"]
# pred.mu_irrign2=raster(pred.mu_irrign2)
# pred.mu_irrign2=mask(pred.mu_irrign2,India_State_Boundary_Bihar_wgs84.prj)
# pred.mu_irrign2_plot=levelplot(pred.mu_irrign2,par.settings=RdBuTheme())
# pred.mu_irrign2_plot

# png("figures/BayesianKriggedMeanValue_Levelplot_irrign2_mu.png")
# pred.mu_irrign2_plot
# dev.off()

# pred.sd_irrign2=pred.grid["pred.sd"]
# pred.sd_irrign2=raster(pred.sd_irrign2)
# pred.sd_irrign2=mask(pred.sd_irrign2,India_State_Boundary_Bihar_wgs84.prj)
# pred.sd_irrign2_plot=levelplot(pred.sd_irrign2,par.settings=RdBuTheme())
# pred.sd_irrign2_plot

# png("figures/BayesianKriggedMeanValue_Levelplot_irrign2_sd.png")
# pred.sd_irrign2_plot
# dev.off()

# writeRaster(pred.mu_irrign2, "figures/mean2.tau.pred.mu.image.tif",overwrite=TRUE) 
# writeRaster(pred.sd_irrign2, "figures/mean2.tau.pred.sd.image.tif",overwrite=TRUE)

# save(cf.mean.sp.pred, file = "cf.mean.sp.pred.RData")
# # predict and probability ------------------------------------------------

# cf.mean.sp.pred.pred.prob=rowSums(cf.mean.sp.pred$p.y.predictive.samples>0)/251
# cf.mean.sp.pred.pred.prob50kg=rowSums(cf.mean.sp.pred$p.y.predictive.samples>0.05)/251
# cf.mean.sp.pred.pred.prob100kg=rowSums(cf.mean.sp.pred$p.y.predictive.samples>0.1)/251
# cf.mean.sp.pred.pred.prob200kg=rowSums(cf.mean.sp.pred$p.y.predictive.samples>0.2)/251
# cf.mean.sp.pred.pred.prob300kg=rowSums(cf.mean.sp.pred$p.y.predictive.samples>0.3)/251
# cf.mean.sp.pred.pred.prob400kg=rowSums(cf.mean.sp.pred$p.y.predictive.samples>0.4)/251


# library(MBA)
# surf <- mba.surf(cbind(coords, tau.hat.cf_irrign), no.X=x.res, no.Y=x.res,extend=TRUE, sp=TRUE)$xyz.est
# #surf <- surf [!is.na(overlay(surf, India_State_Boundary_Bihar_wgs84.prj)),] 
# surf <- as.image.SpatialGridDataFrame(surf) 
# z.lim <- range(surf[["z"]], na.rm=TRUE) 

# pred.grid <- as.data.frame(list(pred.coords,pred.mu=cf.mean.sp.pred.pred.mu,pred.sd=cf.mean.sp.pred.sd,
#                                 pred.prob=cf.mean.sp.pred.pred.prob,pred.prob50kg=cf.mean.sp.pred.pred.prob50kg,
#                                 pred.prob100kg=cf.mean.sp.pred.pred.prob100kg,pred.prob200kg=cf.mean.sp.pred.pred.prob200kg,
#                                 pred.prob300kg=cf.mean.sp.pred.pred.prob300kg,pred.prob400kg=cf.mean.sp.pred.pred.prob400kg))

# coordinates(pred.grid) = c("X", "Y") 
# gridded(pred.grid) <- TRUE 

# writeGDAL(pred.grid["pred.prob"], "figures/mean.tau.pred.prob.image.tif")
# writeGDAL(pred.grid["pred.prob50kg"], "figures/mean.tau.pred.prob50kg.image.tif")
# writeGDAL(pred.grid["pred.prob100kg"], "figures/mean.tau.pred.prob100kg.image.tif")
# writeGDAL(pred.grid["pred.prob200kg"], "figures/mean.tau.pred.prob200kg.image.tif")
# writeGDAL(pred.grid["pred.prob300kg"], "figures/mean.tau.pred.prob300kg.image.tif")
# writeGDAL(pred.grid["pred.prob400kg"], "figures/mean.tau.pred.prob400kg.image.tif")



# library(rasterVis)
# pred.prob100kg=pred.grid["pred.prob100kg"]
# pred.prob100kg=raster(pred.prob100kg)
# levelplot(pred.prob100kg,par.settings=RdBuTheme())

# png("figures/BayesianKriggedMeanValue_plus_probability_100kg_ha_irrign2.png")
# levelplot(pred.prob100kg,par.settings=RdBuTheme())
# dev.off()

6 Certainty equivalent and risk premium

7 References

Antle, J.M. 1983. “Testing the stochastic structure of production: A flexible moment-based approach”. Journal of Business and Economic Statistics 1(3): 192-201. Doi: 10.1080/07350015.1983.10509339.

Antle, J.M. 2010. “Asymmetry, partial moments and production risk.” American Journal of Agricultural Economics 92(5): . Doi: https://doi.org/10.1093/ajae/aaq077.

Di Falco, S., Chavas, J., and Smale, M. 2007. “Farmer management of production risk on degraded lands: the role of wheat variety diversity in the Tigray region, Ethiopia.” Agricultural Economics 36: 147-156. Doi: https://doi.org/10.1111/j.1574-0862.2007.00194.x.

Di Falco, S., and Chavas, J. 2009. “On crop biodiversity, risk exposure, and food security in the highlands of Ethiopia”. American Journal of Agricultural Economics 91(3): 599-611. Doi: https://doi.org/10.1111/j.1467-8276.2009.01265.x.