#######################################################################################
# Charles B. Moss - Drawing Profit Maximizing Output Using the Leontief               #
#######################################################################################
# Heady production function for corn, Heady and Dillon book (1961) p 478              #
#                                                                                     #
#  Input fertilizer price data                                                        #
dta <- read.csv("FertilizerPrices.csv")
#*************************************************************************************#
#  Compute covariance matrix removing the common time trend and interaction with      #
#        phosphorous price                                                            #
ddta <- cbind(dta[,2],dta[,3],dta[,4],dta[,1])
#    mv - mean of relative input price data                                           #
mv   <- mean(ddta[,1:2])
#    vc - variance for complete dataset                                               #
vc   <- var(ddta)
#    vvc - conditional variance matrix                                                #
vvc  <- vc[1:2,1:2]-vc[1:2,3:4]%*%solve(vc[3:4,3:4])%*%vc[3:4,1:2]
#    cvc - Cholesky decomposition for random draw                                     #
cvc  <- chol(vvc)
#*************************************************************************************#
#  Define the Leontief production function                                            #
leon <- function(x) -5.68-0.361*x[1]-0.417*x[2]+6.351*x[1]^0.5+8.5115*x[2]^0.5+
           0.3410*(x[1]^0.5)*(x[2]^0.5)
#  Compute estimated values of the production function given the Leontief structure   #
for (i in 1:nrow(ddta)) {
    #  Setup the profit function for each observation (based on relative prices)      #
    prleon <- function(x) -leon(x) + ddta[i,1]*x[1]+ddta[i,2]*x[2]
    #  Find the level of inputs that maximize profit                                  #
    st     <- optim(c(100,100),prleon)
    #  Form a dataset with maximum profit and individual input levels                 #
    if (i==1) obs <- c(leon(c(st$par[1],st$par[2])),st$par[1],st$par[2]) else
         obs <- rbind(obs,c(leon(c(st$par[1],st$par[2])),st$par[1],st$par[2]))
}
#  Compute the variance of the 'observed' optimums                                    #
obvar <- var(obs)
#  Set the sample variance equal to 0.50 times the diagonal of the 'observed' variance#
rvv   <- chol(0.25*cbind(c(obvar[1,1],0,0),c(0,obvar[2,2],0),c(0,0,obvar[3,3])))
#######################################################################################
#  Randomly draw from the input price and generate profit maximizing output levels    #
for (i in 1:1000) {
    # Randomly draw input prices                                                      #
    pdraw  <- mv + cvc %*% rnorm(2)
    # Setup the profit function for each simulated value - the new objective function #
    prleon <- function(x) -leon(x)+pdraw[1]*x[1]+pdraw[2]*x[2]
    # Choose the level of inputs that maximize profit                                 #
    ss     <- optim(c(5,5),prleon)
    # Randomly draw around the optimal combination of inputs                          #
    #    oo[1] - Random draw around the production function                           #
    #    oo[2] - Random draw around the relative nitrogen price                       #
    #    oo[3] - Random draw around the relative phosphorous price                    #
    oo     <- c(0,ss$par[1],ss$par[2])+rvv%*%rnorm(3)
    # Generate the value on the productin frontier 'leon()' and add the error term    #
    oo[1]  <- leon(cbind(oo[2],oo[3]))+ oo[1]
    # temp observe variance based on deviation from optimal decision                  #
    temp   <-cbind(leon(cbind(ss$par[1],ss$par[2]))-leon(cbind(oo[2],oo[3])),
               ss$par[1]-oo[2],ss$par[2]-oo[3])
    if (i==1) drw <- temp else drw <- rbind(drw,temp)
    # Aggregation into the final dataset                                              #    
    if (i==1) sol <- cbind(oo[1],oo[2],oo[3],pdraw[1],pdraw[2]) else 
      sol <- rbind(sol,cbind(oo[1],oo[2],oo[3],pdraw[1],pdraw[2]))
}
#######################################################################################
#  Setup variables for test regressions                                               #
y <- sol[,1]
n <- sol[,2]
p <- sol[,3]
nn <- sqrt(sol[,2])
pp <- sqrt(sol[,3])
np <- sqrt(sol[,2])*sqrt(sol[,3])
wn <- sol[,4]
wp <- sol[,5]
#*************************************************************************************#
# Leontief regression                                                                 #
lsq <- lm(y~n + p + nn + pp + np)
print(summary(lsq))
# Derived Demand Curve for Nitrogen                                                   #
lsqn <- lm(n~wn + wp)
print(summary(lsqn))
# Derived Demand Curve for Phosphorous                                                #
lsqp <- lm(p~wn + wp)
print(summary(lsqp))
#######################################################################################
# Write Leontief dataset                                                              #
write.csv(sol,file = "Leontief-1.csv")