# Numerical Solvers

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');
```