Recently, I have become interested in numerical analysis withinin macroeconomics. In particular, I have been trying to learn algorithms such as "backwards iteration", the "shooting method" etc, in order to numerically solve an optimal control model, more specifically a growth model.

Could you give me some references for introductory and intermediate numerical analysis textbooks ?

  • 1
    $\begingroup$ I think the "Are they inspiring from other codes or just writing their own code" is unanswerable. The only people that could answer that for you would be the authors of the code. In my opinion, this question would be best without the second and third paragraphs, but that is up to you ultimately. $\endgroup$
    – cc7768
    Aug 1, 2015 at 17:11

3 Answers 3


Read the theory and then try coding. But you're lucky, because I was eager to develope and implement the forward shooting algorithm myself and thus provide it.

Optimal Growth

(See Ben Moll for details.)

Optimal growth model in continuous time reads \begin{align} &\max_c\int^\infty_0 e^{-\rho t}u(c)dt\\ \text{s.t.}~~~& \dot k = f(k) - \delta k - c\\ &c\in[0,f(k)]\\ &k(0) = k_0 \end{align} where \begin{align} f(k)&=Ak^\alpha\\ u(c)&=\frac{c^{1-\sigma}}{1-\sigma}. \end{align} From the first order conditions we derive a system of differential equations \begin{align} \dot k &= f(k) - \delta k - c\\ \dot c &= \frac{c}{\sigma}(f'(k) - \delta - \rho) \end{align}

Fixed points given at \begin{align} \dot k &= 0 \quad \Longrightarrow\quad \tilde k = \left(\frac{\alpha A}{\delta+\rho}\right)^\frac{1}{1-\alpha}\\[2mm] \dot c &= 0 \quad \Longrightarrow\quad \tilde c = f(\tilde k) -\delta \tilde k \end{align}

Since the system is saddlepoint stable there is only one trajectory converging to the unique steady state. Initial stock is given and thus optimal $k^*(0) = k_0$. We aim to solve for $c^*(0)$.

The idea of shooting is to guess $c(0)$ and then iterate over time. If the time paths lead to the steady state we are done, otherwise we have to try a different guess. First of all we need to discretize the differential equations. (Computer works with discrete steps. However, if you like to work with the differential equations try Matlab's boundary value solver bvp4c(). This works out quite easy. Since we have no clue what the solver is doing, we'd like to solve the equations "manually".) The forward difference of some function $\dot x$ is defined as \begin{align} \dot x(t) :\approx \frac{x(t+1)-x(t)}{\Delta t} \end{align}

which gives \begin{align} k(t+1) &= \Delta t(f(k(t)) - \delta k(t) - c(t)) + k(t)\\[2mm] c(t+1) &= \Delta t\frac{c(t)}{\sigma}(f'(k(t)) - \delta - \rho) + c(t). \end{align}

Since the probability of guessing the correct $c^*(0)$ tends to be zero, we need to apply an algorithm which converges towards the true value.


The following code snipped is written in Matlab and presents the main idea.

c_lo = 0;
c_hi = A*k(1)^a;
i = 0;
dist_k = tol + 1;
dist_c = tol + 1;
while (i < I) & (dist_k > tol | dist_c > tol);
c(1) = (c_lo + c_hi)/2;
    for t = 1:T-1;
        k(t+1) = del_t*(A*k(t)^a - d*k(t) - c(t)) + k(t);
        c(t+1) = del_t*c(t)/s*(a*A*k(t)^(a-1) - p - d) + c(t);
        if c(t+1) > A*k(t+1)^a
           c(t+1) = A*k(t+1)^a;
        elseif c(t+1) < 0 
               c(t+1) = 0;
    if k(T) > k_ss & c(T) < c_ss
        c_lo = c(1);
        elseif k(T) < k_ss & c(T) < c_ss
        c_hi = c(1);
        elseif k(T) < k_ss & c(T) > c_ss
        c_hi = c(1);
dist_k = abs(k(T) - k(T-1));
dist_c = abs(c(T) - c(T-1));
i = i + 1;
  • Due to time constraint I can't provide a detaild explanation. I hope it's self-explaining. I may answer specific questions which show some effort.

Boundary Value Problem

As I mentioned above we can easily solve the system by applying Matlab's boundary value problem solver bvp4c. Using $\tilde c = c(T)$ as a terminal condition we have two ODEs and two boundary conditions.

% differential equations with y(1) = k, y(2) = c
dy = @(t,y) [A*y(1)^a - d*y(1) - y(2); ...     % dk/dt;
             y(2)/s*(a*y(1)^(a-1) - p - d)];   % dc/dt
% boundary conditions
bc = @(y0,yT) [y0(1) - k0;     % initial condition
               yT(2) - c_ss];  % terminal condition

solinit = bvpinit(linspace(0,T,10), [k0 c0]);   % initital guess
sol = bvp4c(dy, bc, solinit);                   % call solver
t = linspace(0,T)';  % time axis
y = deval(sol,t)';   % call solution values
k = y(:,1);    % transform variables
c = y(:,2);    % transform variables
  • $\begingroup$ That's really awesome ! Thanks so much for this very very helpful codes and explanations ! $\endgroup$ Sep 3, 2015 at 19:25

You can also find Kenneth Judd's book on numerical analysis here. It is one of the very known and notorious books in that field.


Recursive Macroeconomic Theory isn't a "numerical analysis" book per se, but it has an excellent exposition of the shooting algorithm in a growth model in Chapter 11 of the second edition (I don't think the chapter numbers changed much, but chapter is called "Fiscal Policies in the Nonstochastic Growth Model"). I would recommend starting there.


Ljungqvist, Lars and Sargent, Thomas J. "Recursive Macroeconomic Theory." MIT Press, Cambridge, Mass., 2000.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.