project 3: Unconstrained Optimization Methods
(Due date: Sunday February 25 before 11:30pm in Learn drop box)
This is an individual project. You can help each other, but the final result must be yours.
Question 1
In this question, we want to compute the market demand using a two good market and a one dimen- tional unconstrained Newton method. You have to answer the question using expression and analitical derivatives. The market is composed of 1,000 consumers with the following utility function:
image.png
To generate you 1000 consumers, you will proceed as folows:
set.seed(seed)
s <- runif(1000, 1.5, 5) a1 <- runif(1000, .2, 3) a2 <- runif(1000, 2, 6) I <- runif(1000, 500, 2500) par <- cbind(s=s,a1=a1,a2=a2,I=I)
with the seed number replaced by your student id. Follow the following steps:
a) Write an expression for the utility function and a function that will return either the negative of the utility, the negative of the first derivative or the negative of the second derivarive (with respect to x1) at x1 when prices are p1 and p2, income is I and consumer spends all his income. The function would work as:
getf(ces, list(x1=10, p1=5, p2=7, I=500, a1=1, a2=2, s=3), d=0) ## [1] -200.4395
getf(ces, list(x1=10, p1=5, p2=7, I=500, a1=1, a2=2, s=3), d=1) ## [1] 1.45659
getf(ces, list(x1=10, p1=5, p2=7, I=500, a1=1, a2=2, s=3), d=2)
## [1] 0.04622944
b) Write a function that returns the optimal x1 and x2 using the Newton method. It should solve the following problem:
l <- list(p1=5, p2=7, I=500, a1=1, a2=2, s=3) getsol1(U=ces, x0=3,l=l) ## $sol ## ## 1.734754 70.189461 ## ## $niter ## [1] 7
Here x0 is the starting value for x1 because x2 is replaced by (I −p1x1)/p2, and the optimization is done with respect to x1.
c) The Newton method if not that reliable, expecially when there are restrictions on the possible values for x. For example, if we start too far from the solution, we get the following
getsol1(U=ces, x0=80,l=l,trace=TRUE)
## iteration 1, x1=80, x2=14.29
## iteration 2, x1=-80.73, x2=129.1
## Error in if (crit < tol) break: missing value where TRUE/FALSE needed
We can see that the Newton step is too large, which results in x1 being negative in the second step. That could be a problem because we want to solve 1,000 consumer problems with 1,000 different sets of parameters. Choosing a starting value that will make the algorithm converges for all consumers is impossible. In this step, you will modify your function so that it always converge for any sensible starting value (x0 must be less than I/p1). Here is what you will do. If x10 is the initial x1, the Newton method is:
image.png
where s = fj(x10)/fjj(x10). If x10 is less than the optimal x1, s will be negative and x10 will increase. We don’t want it to increase too much and make x2 negative. If it happens, we want to lower the value of s. If x10 is greater than the optimal x1, s will be positive and x10 will decrease. If it decreases by too much, it will become negative. We can prevent that by decreasing s. Here is the modification you will add to your function: (i) As long as s > x0, then s = s/2. This will prevent x1 to become negative. (ii) As long as s < (x0 − I/p1), then s = s/2. This will prevent x2 to become negative. Something like the following (before computing the new x1) will do the job:
while (s>x0 | s<(x0-l$I/l$p1)) s <- s/2
The new function should now work for any x0:
getsol2(U=ces, x0=80,l=l)
## $sol
## x1 x2
## 1.734754 70.189461
##
## $niter ## [1] 10
getsol2(U=ces, x0=99,l=l)
## $sol
## x1 x2
## 1.734754 70.189461
##
## $niter ## [1] 9
Once your function is stable, remove the argument x0 and set the starting value to a sensible value inside the function. Remember that x2 must be positive.
getsol(U=ces, l)
## $sol
## x1 x2 ## 1.734754 70.189461 ##
## $niter ## [1] 8
c) Write a function that computes the total quantity demanded by your market (the sum of all 1,000 consumers) for a given set of prices p = p1, p2 j. Try to make it work for any number of consumers. You would get something like:
library(parallel)
marketD(p=c(3,6), par, U=ces, tol=1e-7, maxit=100, mc.cores=4)
## Market demand
## *********************** ## Number of consumers: 1000 ## p1 = 3, p2 = 6
## x1 = 48968.99, x2 = 222166.76
## ***********************
Of course, you should all get different values as your par will be different. You are free to create a particular object with its print method if you want, but it is not required.
d)
Write a function that plots the market demand of either x1 for a given p2 or x2 for a given p1. You would get something like:
image.png
Hint: If you want to smooth a function with a few points, you can use spline(). For example:
image.png
e) Compare the effect of a 20% income on both demands when the other price is 5. I am expecting something like
image.png
Question 2
Estimate the ρ from the second project using the Newton method and numerical derivative. To do so, you need to have a function that depends on ρ, and returns the SSR. You can use part of your codes from the second project. Basically, for each ρ, you need to estimate the model:
image.png
set.seed(112233)
u <- arima.sim(n=200, model=list(ar=.7)) x1 <- rnorm(200) x2 <- rnorm(200) y <- 1+x1-x2+u X <- cbind(Intercept=1,x1,x2) (ssr <- SSR(rho=.5, y, X, "CO")) ## [1] 202.0671 (ssr <- SSR(rho=.5, y, X, "PW")) ## [1] 202.0674
After, you need a function that computes the first and second derivative of SSR(ρ).
dssr <- dSSR(rho=.5, y, X, "CO",h=.001)) ## [1] -129.2444 (ddssr <- ddSSR(rho=.5, y, X, "PW",h=.001)) ## [1] 709.8033
Then, write the algorithm to minimize SSR. You should get:
getrho(rho0=.1, y, X, "CO")
## $rho
## [1] 0.6819157 ##
## $SSR
## [1] 190.3081 ##
## $niter ## [1] 3
getrho(rho0=.1, y, X, "PW") ## $rho
## [1] 0.6819152
##
## $SSR
## [1] 190.3086 ##
## $niter ## [1] 3
