Optimum.R

#################################

#        Optimization           #

# 1) function f maps the possible
#  variable domain to the objective
#  function u(var1, var2, var3)

# 2) it compute the value of u and 
#  return the highest value of u 
#  and the optimizing variable set.

# 3) Compute the value function f
#  given the full range of indep
#  -endent variables 'w' and 'b'.

# 4) plot the surface of the value
#  function f using 'w' and 'b'.

# 5) identify income thresholds
#  for corner solutions in u().

# 6) measure the margional u of
#  health ?? at different wealth
#  and the WTP for health. 

#################################



# PART 1
# FUNCTION f


f <- function(w, h){
        n <- 100
        m <- NULL
        u <- function(var1, var2, var3){
                (var1)^(1/2) + ((var2)^(1/2) )
                * (var3+h)/(var3+h+1)
        }
        
        x <- seq(0, w, w/n)
        v <- w - x
        
        for (i in c(1:(n+1))){
                        y <- seq (0, v[i], w/n)
                        z <- v[i] - y
                        s <- cbind(x[i],y,z)
                        m <- rbind(m,s)
                } 
        
        uu <- NULL   
        
        if (w != 0) {
                m <- as.matrix(m)
                for (i in 1:nrow(m)) { 
                        
                um <- data.frame()
                uu <- c(uu, u(m[i, 1], m[i, 2], m[i,3]))
                        
                if (i == nrow(m)) {
                        um <- data.frame("utility" = uu, 
                                        "cons" = m[, 1], 
                                        "save" = m[, 2],
                                        "health" = m[, 3])
                }
                
                }} else { 
                        um <- data.frame("utility" = 0, 
                                         "cons" = 0, 
                                         "save" = 0, 
                                         "health" = 0)
                }
        
        solution <- subset(um, utility == max(um$utility))
        return(solution[1])
        
}


# PART 2
# PLOT FUNCTION f(w, b)

# Generate full range of wealth and health 
# and map them into the optimal utility f()
w = NULL
b = NULL 
w = seq(1, 150, 1)
b = seq(1, 10, 9/(length(w)-1))

# Given the range of wealth w and 
# initial health level b, we can 
# compute the optimal surface of 
# u(c, s, h) in f(w, b). 
# we save the output in string z
# Assume the range for w is (0, 150)
# range for b is (0, 10)

z = NULL
for (j in 0:150) {
        for (i in 0:150){
                z <- c(z, as.numeric(f(w[i], h[j])))
        }
}

z.mat = matrix(z, 150, 150) # this is the output matrix

# PART 2.1
# the overall value surface looks fine...
persp(w,b,z.mat, phi=45,theta=-45,col="yellow",
      shade=.00000001,ticktype="detailed")


# PART 2.2
# to make the visuals for marginal analysis
# scale of wealth (w) and initial health (b) 
# needs to be more comparable.

par(mfrow = c(1,2))
persp(w[1:10],b[1:100],z.mat[1:10, 1:100], 
      phi=10,theta=-45,col="yellow", shade=.00000001,
      ticktype="detailed")
persp(w[1:100],b[1:100],z.mat[1:100, 1:100], 
      phi=10,theta=-45,col="yellow", shade=.00000001,
      ticktype="detailed")

# to make the margins more visible

png(filename="Health and Wealth Effects.png", height=400, width=900, 
    bg="white")

par(mfrow =c(1, 3))
        plot(w, z.mat[,1], ylim = range(z.mat), 
             type ='l', col = "red", lwd = 2, lty = c(2:1),
             xlab = "Wealth", ylab = "Utility", cex = 2,
             main = "Effect of Health on Expected Utility")
        lines(w, z.mat[,150], ylim = range(z.mat), 
             type ='l', col = "forestgreen", lwd = 2
             )
        legend("topleft", legend = c("high", "low"), 
               title = "Health level", col= c("forestgreen", "red"), 
               lty=1:2, lwd=2, bty="n")
        plot(h, z.mat[1,], 
             type ='l', col = "purple", lwd = 2, cex = 2,
             xlab = "Health", ylab = "Utility", ylim = c(1.0, 1.6),
             main = "Wealth Effect on the Utility of Health")
        legend("topleft", legend = "Low", cex=0.8, 
               title = "Wealth level", col= "purple", 
               lty=1, lwd=2, bty="n")
        plot(h, z.mat[150,], ylim = c((16.0), (16.6)), 
             type ='l', col = "blue", lwd = 2, cex = 2,
             xlab = "Health", ylab = "Utility",
             main = "Wealth Effect on the Utility of Health"
             )
        legend("topleft", legend = "High", cex=0.8, 
               title = "Wealth level", col= "blue", 
               lty=1, lwd=2, bty="n")
       dev.off()