Numerical Solvers

This section reviews numerical solvers in MATLAB.

A very useful feature of Matlab is its ability to find solutions to nonlinear equations and, more generally, systems of nonlinear equations.

Finding the root of a nonlinear function

The function fzero searches for a point $x$ where target_function$(x)=0$. It stops searching once it finds a point where the function changes sign. This can be highly useful for example to find solutions to first order conditions.

The syntax for using fzero is

[x,fval,exitflag] = fzero(@function,x0,options)

where $x$ is the solution, fval the value of the function at the solution, exitflag gives the exit condition, function is our function of interest, $x0$ is the starting value (can be a single point or an interval), and options specifies the options of fzero to be used.

Say, for example, if we want to find the root of the sine function close to 3 we write

options = optimset('fzero');
options = optimset(options,'Display','iter');
[x,fval,exitflag] = fzero(@sin,3,options)

Likewise, instead of looking close to 3 we can specify a range by writing

[x,fval,exitflag] = fzero(@sin,[2.5,3.5],options)

However, if the function of interest has hyperparameters we need to make an extra step. Take our univariate function from the previous section, $f(x) = a + bx+cx^2$. The first order conditions are then $b+2cx = 0$ which we can solve with fzero but we need to specify the values of $b$ and $c$. We can write

options = optimset('fzero');
options = optimset(options,'Display','iter');
myfun = @(x,b,c) b+2*c*x;
b = 2;
c = 1.5;
fun = @(x) myfun(x,b,c)
[x,fval,exitval] = fzero(fun,-1,options)

A short cut could be to write

options = optimset('fzero');
options = optimset(options,'Display','iter');
b = 2;
c = 1.5;
[x,fval,exitval] = fzero(@(x) b+2*c*x,-1,options)

Solving systems of nonlinear equations

Often the problems we face are not simply a function of a single variable but are rather systems of nonlinear equations. To solve these problems we can use the fsolve function which solves $F(x)=0$ for $x$ (which is a vector or matrix), where $F(x)$ is a function returning a vector value. The syntax is comparable to that of fzero, that is

[x,fval,exitflag] = fsolve(@function,x0,options)

where $x$ is the solution, fval the value of the function at the solution, exitflag gives the exit condition, function is our function of interest, x0 is the starting value, and options specifies the options of fsolve to be used.

Say that we have a system of two equations in two variables:

$exp(-x) = y^2$

$y cos(x) = 1$

To solve this system we first need to convert the system to the form $F(X)=0$ and write a function that calculates the left-hand side of the the equation:

function F = d2roots(z)
F(1) = exp(-z(1)) - z(2)^2;
F(2) = z(2)*cos(z(1)) - 1;
end

Then we can pass the function to fsolve to obtain the solution:

options = optimset('fsolve');
options = optimset(options,'Display','iter');
z0 = [0,0];
[z,fval,exitval] = fsolve(@d2roots,z0,options)

Solving a New Keynesian model using fsolve

The basic New Keynesian model consists of three core equations:

$\pi_t = \beta E_t \pi_{t+1}+\kappa y_t$ (The NK Phillips curve)

$y_t = E_t y_{t+1} - \sigma^{-1}(i_t - E_t \pi_{t+1}-\rho+z_t)$ (The Dynamic IS equation)

$i_t = \rho + \phi_{\pi}\pi_t + \phi_y y_t$ (Monetary policy)

where $\pi_t$ denotes the inflation rate, $y_t$ the output gap, $i_t$ the nominal interest rate and $z_t$ is a "financial shock" which evolves according to the process $z_t = \rho_z z_{t-1} + \epsilon_t$ where $\epsilon_t$ is an iid random variable with zero mean and standard deviation $\sigma_{\epsilon}$.

We consider a finite horizon version of the model and use fsolve to solve the system and give us the impulse response of a one standard deviation decrease in $\epsilon_t$ (i.e., $\epsilon_0 = -\sigma_{\epsilon}$, $\epsilon_t = 0 \ \ \forall t>0$) for the parameter values

$\beta$	$\kappa$	$\sigma$	$\phi_{\pi}$	$\phi_y$	$\rho_z$	$\sigma_{\epsilon}$	$\rho$
0.99	0.17	1	1.5	1/8	0.8	0.02	$-log(\beta)$

We begin by setting up our target function (recall that our target function needs to have all the choice parameters as the first input):

function out = BNK_mdl(x,T,nvar,zshock)
        % Parameters
        betta = 0.99;        % discount factor
        rho = -log(betta);   % steady state interest rate
        kappa = 0.17;       % slope of the NK Phillips curve
        sigma = 1;          % CRRA coefficient
        phipi = 1.5;        % CB weight on inflation
        phiy = 1/8;         % CB weight on output

        % Define variables
        x = reshape(x,T,nvar);
        pi = x(:,1);
        y = x(:,2);
        int = x(:,3);

        % Create the forward terms
        pi_next = [pi(2:end);0];    % The steady state pi=0 is the terminal condition
        y_next = [y(2:end);0];      % The steady state y=0 is the terminal condition

        % Equations being solved
        z(:,1) = -pi + betta*pi_next + kappa*y;                          % NK Phillips curve
        z(:,2) = -y + y_next - 1/sigma.*(int - pi_next - rho + zshock); % DIS equation
        z(:,3) = -int + rho + phipi*pi + phiy*y;                        % Monetary Policy Rule

        out=z(:); % stack the columns to output a column vector
end

And now, to solve the model and to plot the impulse responses we run the following:

%% Solve the model
% Set the options for the root finding algorithms
options = optimset('fsolve');
option=optimset(options,'MaxFunEvals',200000,'MaxIter',10000,'TolFun',1e-8,'TolX',1e-8, 'Display','iter');

% Specify some variables
T = 30; % Number of periods
nvar = 3; % Number of variables
varnames = char('Inflation','Output Gap','Int.rate');

% Create the path of the financial shock
rhoz = 0.8;         % persistence of financial shock
sige = 0.02;        % s.d. of innovation of financial shock
zshock = NaN(T,1);
zshock(1,1) = -sige;
for i = 2:T
    zshock(i,1) = rhoz*zshock(i-1,1);
end

%% Solve model using
xstart = zeros(nvar*T,1);                   % initial guess
solpath = fsolve(@(x) BNK_mdl(x,T,nvar,zshock),xstart(:),option); % find solution

%% Plot solutions
solpath = reshape(solpath,T,nvar);
% Plot solutions using fsolve
figure('Name','Using fsolve')
for i=1:nvar
    subplot(2,2,i)
    plot([0:T-1], solpath(:,i),'.-');
    title(varnames(i,:));
    xlabel('Period');
    ylabel('Impact');
end
subplot(2,2,4)
plot([0:T-1], zshock,'.-');
title('Financial shock');
xlabel('Period');
ylabel('Impact');