# Author: Mattias Villani, Sveriges Riksbank and Stockholm University. # E-mail: mattias.villani@riksbank.se # Script to illustrate numerical maximization of the function # f(x,y)=(x-a)^2+a*y^2-2*x*sin(y). Note that the above expression is really # a function of x,y and a, but I write as f(x,y) since I will maximize # with respect to x and y for a given a. # First define the function f. Note that the minimization routine in R # forces us to group all the variables which we intended to maximize # with respect to in a vector z=(x,y), so that x=z[1] and y=z[2]. f<-function(z,a){ t=(z[1]-a)^2+a*z[2]^2-2*z[1]*sin(z[2]) return(t) } # Initial values of z Init=c(5,-4) # Calling the optimization routine. Note the argument a=3 OptimResults<-optim(Init,f,hessian=TRUE,a=3) XOpt=OptimResults\$par[1] # x-value at the optimum YOpt=OptimResults\$par[2] # y-value at the optimum # Evaluate the function over two grids (one for x and one for y) # to see if we indeed have a maximum XGrid=seq(2,5,0.1) YGrid=seq(0.5,1.5,0.1) FunctionXGrid=matrix(NA,length(XGrid)) FunctionYGrid=matrix(NA,length(YGrid)) a=3 XCount=0; # Evaluating the function f(x,y) over The x-grid, with y kept fixed at YOpt for (x in XGrid){ XCount=XCount+1; FunctionXGrid[XCount]=f(c(x,YOpt),a) } YCount=0; # Evaluating the function f(x,y) over The y-grid, with x kept fixed at XOpt for (y in YGrid){ YCount=YCount+1; FunctionYGrid[YCount]=f(c(XOpt,y),a) } # Plotting the perturbated function par(mfrow=c(1,2)) # Opens a graphic windows consisting of two subplots which # are filled row by row plot(XGrid,FunctionXGrid,type='l') # The extra argument type='l' makes it a line plot title('Function over the x-grid for y=YOpt') plot(YGrid,FunctionYGrid,type='l') title('Function over the y-grid for x=XOpt')