Introduction to MATLAB

Session 5: Introduction to Numerical Optimization and MATLAB

Solvers

Danila Smirnov

Barcelona School of Economics

September 2022
Introduction
We want to find the solution of a system of equations or the max/min of
a function
ˆ Find the roots or zeroes of a function: f (x ∗ ) = 0
ˆ Find the maximize/minimize a function: f 0 (x ∗ ) = 0

Examples
ˆ Find the equilibrium price of an economy
ˆ Maximize a log-likelihood function

Numerical methods will help us to find a numerical solution for the

problem at hand

Key
ˆ Because we are numerically solving a problem, we have to accept
some tolerance for the solution
ˆ We will rely on iterative algorithms: the starting point matters
ˆ Local solutions
2/23
Example
Where are the zeroes? Where is the function maximized/minimized?

f (x) = x cos(x 2 )

-1

-2

-3
-1 -0.5 0 0.5 1 1.5 2 2.5 3

3/23
Rootfinding problems
Bisection

ˆ Algorithm to find the root of a continuous real-valued function

ˆ Based on the Intermediate Value Theorem: if f is continuous, and
f (a) and f (b) have different signs, then f must have at least one
root x in [a, b]
ˆ Pros: very robust algorithm
ˆ Cons: only applicable to one-dimensional problems

4/23
Rootfinding problems
Bisection: Graphical example

Source: https://orionquest.github.io/Numacom/bisection.html

5/23
Rootfinding problems I
Bisection

The pseudo-code for bisection is as follows

1. Determine the function f and the interval in which you want to find
the zero (say [a, b]). Important! check that the sign of f (a) is
different to the sign of f (b)
◦ Save them as s_min and s_max
2. Set an initial guess for the solution inside [a, b]. For example,
x 0 = (a+b)
2
3. Set a tolerance and an initial error term (corresponding to the length
of the interval)
4. Initiate the iterative step
4.1 Evaluate the sign of the function at the candidate solution f (x 0 )
If sign f (x 0 ) = s_min, increase the initial guess
Else, decrease the initial guess
4.2 Iterate until the length of the interval is smaller than the tolerance
go to MATLAB now! (example 1)

6/23
Rootfinding problems
Newton’s Methods

Iterative method based on successive linearization

f (x) ≈ f (x k ) + f 0 (x k )(x − x k ) = 0 (1)

Basic algorithm
1. Guess x (0) for the root of f . The super-script 0 denotes the iteration
number
2. Using (1), update the guess with

x (k+1) = x (k) − [f 0 (x (k) )]−1 f (x (k) ) (2)

3. Iterate until x (k+1) and x (k) are “close”

Any “problem” with this algorithm?

7/23
Rootfinding problems
Quasi-Newton’s Methods

The previous algorithm requires the Jacobian (f 0 ) of the function

Quasi-Newton methods replace the Jacobian by variants of a numerical

derivative

f (x (k) ) − f (x (k−1) )
f 0 (x (k) ) ≈
x (k) − x (k−1)

The iteration rule now reads as

x (k) − x (k−1)
x (k+1) = x (k) − f (x (k) ) (3)
f (x (k) ) − f (x (k−1) )

Any “problem” with this algorithm?

8/23
Rootfinding problems I
MATLAB functions

Root of non-linear function

[x,fval,exitflag,output] = fzero(‘fun’,x0,options)

Root of a system of non-linear equations

[x,fval,exitflag,output] = fsolve(‘fun’,x0,options)

Input
ˆ ‘fun’: function we want to find a zero
ˆ x0: initial guess for the solution.
ˆ Depending on the problem, you can set different options
Output
ˆ x: zero of equation ‘fun’
ˆ fval: value of the function at the optimal point
ˆ exitflag: integer encoding the exit condition
ˆ output: information about the root-finding process

9/23
Rootfinding problems II
MATLAB functions

Important!
When using fzero, you a single equation and must provide a single
initial starting point. When using fsolve, you have a system of N
equations and have to provide a N-dimensional vector as starting point.

10/23
Exercise 1

Write a function f1 defining f (x) = x cos(x 2 )

1. Find the zero of this function using fzero and initial guess x 0 = 1
2. Repeat using x 0 = 2 as initial guess
3. Repeat using fsolve
4. How does the solution changes with the algorithm used? What
about the initial guess?

11/23
Rootfinding problems
Setting options

For fzero use optimset:

ˆ Display: level of display. For example ’off’ displays no output,
while ’iter’ displays output at each iteration
ˆ TolX: termination tolerance on x.
Example: options = optimset(’Display’,’off’,’TolX’,1e-8)

For fsolve use optimoptions and assign to fsolve

ˆ Display: same as before
ˆ MaxIter: maximum number of iterations allowed (positive integer)
ˆ TolFun: termination tolerance on the function value (positive scalar)
Example:
options = optimoptions(’fsolve’,’Display’,’off’,’TolFun’,1e-8,’MaxIter’,1000)

Note! you can set more options. See help fzero and help fsolve for
details

12/23
Rootfinding problems
MATLAB functions: Passing arguments

Suppose you have a function

function y = myfunction(x,a,b)
Where a and b are arguments/parameters of the function.

You can use fzero/fsolve by calling myfunction as an anonymous

function (recall Session 3)
fzero(@(x) myfunction(x,a,b),x0,options)

13/23
Numerical optimization
Two commonly used methods

Grid search
ˆ Generate a fine grid of points and evaluate the function at each one.
Choose x ∗ as the value that generates the highest/lowest f (x ∗ )

Derivative based
ˆ At the optimum, the first derivative of the function is zero. As in
root-finding methods, this is an iterative procedure that uses
numerical derivatives

14/23
Numerical optimization
Grid search

Basic steps
1. Generate a first rough grid with n points and evaluate the function
2. Select the point that maximizes/minimizes the function on this grid,
x (k)
3. Refine the initial grid around the optimal point x (k)
4. Select the new point that optimizes the function
5. Repeat steps 3-4 until the optimized values of consecutive iterations
are “close”

Even though this method is clear and intuitive, it computationally

expensive

go to MATLAB now! (example 2)

15/23
Numerical optimization
Newton-Raphson

Same spirit as for rootfinding algorithm: linear approximation

g (x) ≈ g (x k ) + g 0 (x k )(x − x k )

Solving the first-order condition, g (x) = f 0 (x)

f 0 (x k ) + f 00 (x k )(x − x k ) = 0

Updating rule

x k+1 = x k − [f 00 (x k )]−1 f 0 (x k )

16/23
Numerical optimization
MATLAB functions

Minimum of unconstrained multivariable function using derivative-free

method
[x,fval,exitflag,output] = fimnsearch(fun,x0,options)

Minimum of unconstrained optimization

[x,fval,exitflag,output] = fminunc(fun,x0,options)

Input/output and passing additional arguments is similar to fzero and

fsolve

Note! for fminsearch you can set similar options as with fzero (using
optimset). See help fminsearch.

Important!
When using fminsearch or fminunc you can have a N-dimensional
vector as starting point, but the output of function fun must be a scalar

17/23
Exercise 2

Recall the function f (x) = x cos(x 2 )

1. Find the maximum of this function using fminsearch. Use as initial
guess x 0 = 1
2. Repeat the previous step but using as initial guess x 0 = 3. Any
difference?
3. Find the maximum of this function using fminunc. Use as initial
guess x 0 = 1

18/23
Numerical optimization
Tips

Suppose you want to obtain θ which must satisfy some constraint

Get an unconstrained value ψ and transform as follows

Example Constraint Transformation

Variance θ>0 θ = ψ2

θ = exp(ψ)
1
Probability θ ∈ (0, 1) θ= 1+exp(ψ −1 )
ψ
Stationary autoregressive parameter θ ∈ (−1, 1) θ= 1+|ψ|
1
Shares θ1 , θ2 ≥ 0 θ1 = 1+exp(ψ)
exp(ψ)
θ1 + θ2 = 1 θ2 = 1+exp(ψ)

19/23
Exercise 3
Consider a Cournot duopoly model, in which the inverse demand is:
P(q) = q 1/η

Each firm i = 1, 2 have costs: C (qi ) = 12 ci qi2

Profits of firm i are: πi (q1 , q2 ) = P(q1 + q2 )qi − Ci (qi )

The first order condition reads as

∂πi
= P(q1 + q2 ) + P 0 (q1 + q2 )qi − Ci0 (qi ) = 0
∂qi
Thus the market equilibrium outputs, q1 and q2 , are the roots of the two
nonlinear equations

fi (q) = (q1 + q2 )−1/η − (1/η)(q1 + q2 )−1/η−1 qi − ci qi = 0, i = 1, 2

20/23
Exercise 3

Assume η = 1.6, c1 = 0.6 and c2 = 0.8 and find the equilibrium of this
economy
ˆ Create a function called cournot_res that receives as input
1. The initial guess for the vector of equilibrium output
2. The values of parameters (η, c1 , c2 )
and takes as output the “residual” of the equilibrium conditions
ˆ Use one the solvers covered previously (which one? why?)

How the equilibrium output of each firm changes with the elasticity η?
ˆ Set a linear space of 100 values between 0.6 and 5 for η
ˆ Solve the model for each value on the grid
ˆ Plot the equilibrium output of each firm against the elasticity

21/23
Exercise 4
Suppose you want to estimate the following model

yi = α + βxi + ui

ˆ xi : regressor of interest
ˆ yi : dependent variable
ˆ ui ∼ N(0, σ 2 ): error term, independent of xi and independent across
observations

The log-likelihood reads as

n
n 1 X
` = − log(2π) − n log(s) − 2 (yi − α − βxi )2 (4)
2 2s
i=1

s 2 : estimator of the variance error

22/23
Exercise 4

Load the dataset us_data.xlsx

ˆ spread: difference between the 10-year yield and the 2-year yield
ˆ growth: growth rate of GDP of the US

Estimate the parameters (α, β, s 2 ) by Maximum Likelihood Estimation

(MLE)
ˆ Create a function us_mle that receives as input
1. The initial guess for the vector of parameters
2. The regressor and the independent variable
and takes as output the log-likelihood (4)
Hint! Do you have to impose any constraint in the estimation?

Compare the solution with Ordinary Least Squares: βb = (X 0 X )−1 X 0 Y

23/23