<p>You can also find Kenneth Judd's book on numerical analysis <a href="https://mitpress.mit.edu/books/numerical-methods-economics" rel="nofollow">here</a>. It is one of the very known and notorious books in that field.</p>

CC BY-SA 3.0
2015-08-23T21:03:56.363Z
2015-08-23T21:03:56.363Z
user1905
<p>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.</p>

<p><strong>References</strong></p>

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

CC BY-SA 3.0
2015-08-01T17:09:40.197Z
2015-08-01T17:09:40.197Z
cc7768
<p>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.</p>

<h2>Optimal Growth</h2>

<p>(See <a href="http://www.princeton.edu/~moll/ECO503Web/Lecture3_ECO503.pdf" rel="nofollow">Ben Moll</a> for details.)</p>

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

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

<p>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)$. </p>

<p>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 <code>bvp4c()</code>. 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}</p>

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

<p>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. </p>

<h2>Shooting</h2>

<p>The following code snipped is written in Matlab and presents the main idea. </p>

<pre><code>c_lo = 0;
c_hi = A*k(1)^a;
i = 0;
dist_k = tol + 1;
dist_c = tol + 1;
while (i &lt; I) &amp; (dist_k &gt; tol | dist_c &gt; 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) &gt; A*k(t+1)^a
           c(t+1) = A*k(t+1)^a;
        elseif c(t+1) &lt; 0 
               c(t+1) = 0;
        end
    end
    if k(T) &gt; k_ss &amp; c(T) &lt; c_ss
        c_lo = c(1);
        elseif k(T) &lt; k_ss &amp; c(T) &lt; c_ss
        c_hi = c(1);
        elseif k(T) &lt; k_ss &amp; c(T) &gt; c_ss
        c_hi = c(1);
    end
dist_k = abs(k(T) - k(T-1));
dist_c = abs(c(T) - c(T-1));
i = i + 1;
end
</code></pre>

<ul>
<li>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.  </li>
</ul>

<h2>Boundary Value Problem</h2>

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

<pre><code>% 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
</code></pre>

CC BY-SA 3.0
2015-09-03T19:25:47.133Z
That's really awesome ! Thanks so much for this very very helpful codes and explanations !
optimal control
CC BY-SA 3.0
2015-09-03T13:54:23.603Z
2015-09-03T17:13:19.960Z
2015-09-03T17:13:19.960Z
clueless
<p>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. </p>

<p>Could you give me some references for introductory and intermediate numerical analysis textbooks ? </p>

CC BY-SA 3.0
2015-08-01T17:11:04.740Z
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.
cc7768
<macroeconomics><numerical-methods>
Numerical analysis in economics
optimal control
# Numerical analysis in economics

(Asked by: username_0 on 2015-08-01)

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 ?

**Tags:** macroeconomics, numerical-methods
> **Comments:**
> **username_1**: 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.  

---

## 3 answers

### Answer by username_3 on 2015-09-03. Score: 2 | Accepted answer

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. 
Shooting
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;
        end
    end
    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);
    end
dist_k = abs(k(T) - k(T-1));
dist_c = abs(c(T) - c(T-1));
i = i + 1;
end


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

> **Comments:**
> **username_0**: That's really awesome ! Thanks so much for this very very helpful codes and explanations !  

---

### Answer by username_2 on 2015-08-23. Score: 3

You can also find Kenneth Judd's book on numerical analysis [here](https://mitpress.mit.edu/books/numerical-methods-economics). It is one of the very known and notorious books in that field.


---

### Answer by username_1 on 2015-08-01. Score: 1

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.

**References**

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


---

