
\documentclass[12pt,a4paper,titlepage]{article}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\usepackage{amsmath}
\usepackage{amsmath}
\usepackage{amsmath}
\usepackage{amsmath}
\usepackage{amsmath}
\usepackage{amsmath}
\usepackage{amsmath}
\usepackage{amsmath}
\usepackage{amsmath}
\usepackage{amsmath}
\usepackage{amsmath}
\usepackage{amsmath}
\usepackage{amsmath}
\usepackage{amsmath}
\usepackage{amsmath}
\usepackage{amsmath}
\usepackage{amsmath}
\usepackage{amsmath}
\usepackage{amsmath}
\usepackage{amsmath}
\usepackage{amsmath}
\usepackage{amsmath}
\usepackage{amsmath}
\usepackage{amsmath}
\usepackage{amsmath}
\usepackage{amsmath}
\usepackage{amsmath}
\usepackage{amsmath}
\usepackage{amsmath}
\usepackage{amsmath}
\usepackage{graphicx}
\usepackage{amsmath}
\usepackage{harvard}

\setcounter{MaxMatrixCols}{10}
%TCIDATA{OutputFilter=LATEX.DLL}
%TCIDATA{Version=5.50.0.2960}
%TCIDATA{<META NAME="SaveForMode" CONTENT="1">}
%TCIDATA{BibliographyScheme=Manual}
%TCIDATA{Created=Mon Nov 05 11:50:44 2001}
%TCIDATA{LastRevised=Tuesday, December 14, 2010 09:51:23}
%TCIDATA{<META NAME="GraphicsSave" CONTENT="32">}
%TCIDATA{<META NAME="DocumentShell" CONTENT="Journal Articles\Standard LaTeX Article (Harvard)">}
%TCIDATA{Language=American English}
%TCIDATA{CSTFile=LaTeX article (bright).cst}
%TCIDATA{ComputeDefs=
%$\delta =.2$
%$h=2$
%$q=1/5$
%1$s_{e}=1$
%1$s_{u}=-2$
%$\alpha =1/2$
%}


\newtheorem{theorem}{Theorem}
\newtheorem{acknowledgement}[theorem]{Acknowledgement}
\newtheorem{algorithm}[theorem]{Algorithm}
\newtheorem{axiom}[theorem]{Axiom}
\newtheorem{case}[theorem]{Case}
\newtheorem{claim}[theorem]{Claim}
\newtheorem{conclusion}[theorem]{Conclusion}
\newtheorem{condition}[theorem]{Condition}
\newtheorem{conjecture}[theorem]{Conjecture}
\newtheorem{corollary}[theorem]{Corollary}
\newtheorem{criterion}[theorem]{Criterion}
\newtheorem{definition}{Definition}
\newtheorem{example}[theorem]{Example}
\newtheorem{exercise}[theorem]{Exercise}
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{notation}[theorem]{Notation}
\newtheorem{problem}[theorem]{Problem}
\newtheorem{proposition}[theorem]{Proposition}
\newtheorem{remark}[theorem]{Remark}
\newtheorem{result}[theorem]{Result}
\newtheorem{solution}[theorem]{Solution}
\newtheorem{summary}[theorem]{Summary}
\newenvironment{proof}[1][Proof]{\textbf{#1.} }{\ \rule{0.5em}{0.5em}}
\input{tcilatex}
\begin{document}

\title{MathII at SU and SSE}
\author{John Hassler \\
%EndAName
IIES, Stockholm University}
\maketitle

\section{Introduction\protect\footnote{%
I am grateful for comments and corrections provided by Dirk Niepelt.}}

This course is about dynamic systems, i.e., systems that evolve over time.
The analysis of the dynamic evolvement of economic systems is of core
importance in many areas of economics. Growth, business cycles, asset
pricing, and dynamic game theory are just a few examples.

\subsection{Solving a simple dynamic system}

Very often, our economic models provide a difference or differential
equation for the endogenous variables. Take a very simple example; an
arbitrage asset pricing model. Suppose there is a safe asset, a bond, that
provides a return $r$ every period. Suppose a share, giving rights to
dividend flow of $d$ per period$,$ is introduced. Now, arbitrage theory says
that the share should also yield a return $r$ in equilibrium. Defining the
price on the share as $p$, arbitrage theory thus implies,

\begin{equation}
\frac{p_{t+1}+d}{p_{t}}=1+r.
\end{equation}

This is a simple difference equation, 
\begin{equation}
p_{t+1}=\left( 1+r\right) p_{t}-d.
\end{equation}

Note that we $p_{t}$ is an \emph{endogeneous} variable in this system while $%
r$ and $d$ are \emph{exogeneous}. This would be the case also if $d$ and $r$
varies over time, cases we will study below.

One straightforward way of solving it is to substitute forward or backward,
e.g., noting that 
\begin{eqnarray}
p_{t} &=&\left( 1+r\right) p_{t-1}-d \\
&=&\left( 1+r\right) \left( \left( 1+r\right) p_{t-2}-d\right) -d \\
&=&\left( 1+r\right) ^{2}p_{t-2}-d\left( 1+\left( 1+r\right) \right)
\end{eqnarray}%
and so on. A more general approach is to first characterize all possible
paths consistent with the law-of-motion. Here, this is quite simple. You
will learn that set of possible paths is 
\begin{equation}
p_{t}=c\left( 1+r\right) ^{t}+\frac{d}{r}
\end{equation}%
for any constant $c.$ As we see, there is an infinite number of solutions,
i.e., we need more information. If, for example, we know that the value of $%
p_{0}$, we can solve for the constant 
\begin{eqnarray}
c &=&p_{0}-\frac{d}{r} \\
&\rightarrow &p_{t}=\left( p_{0}-\frac{d}{r}\right) \left( 1+r\right) ^{t}+%
\frac{d}{r}.
\end{eqnarray}

In finance, the solution $p=\frac{d}{r}$ is called the fundamental solution,
and we see that if $r>0,$ the solution explodes as $t$ goes to infinity if $%
p_{0}\neq \frac{d}{r}.$ This gives us another way of solving the difference
equation. Suppose, we have reason to believe that the solution should remain
bounded. Then, if $r>0$, the only solution left is when $c=0.$ Note that $%
\left( 1+r\right) $ is called the root of the system. Note, the importance
of whether the root is bigger or smaller than unity (in absolute values).

We will also work in continuous time. Then,\ I usually use put the time
variable in parenthesis and use the dot symbol to indicate time derivatives.
In continuous time, non-existence of arbitrage means that capital gains,
i.e., the change in the price per unit of time plus dividends per unit of
time should equal the opportunity cost, i.e., interest rate on the price of
the assets. Thus, non-existence of arbitrage implies 
\begin{equation}
\dot{p}\left( t\right) +d=rp\left( t\right) .
\end{equation}%
This is a simple linear first-order differential equation. The set if
solution is:%
\begin{equation}
p\left( t\right) =ce^{rt}+\frac{d}{r}.
\end{equation}

Also here, we have a term that is explosive if $r>0.$

Later in the course, we will learn how to solve more complicated dynamic
systems, involving, e.g., several endogenous variables and varying
parameters.

\subsection{Two approaches to Dynamic Optimization}

The second part of the course, is to solve maximization problems in dynamic
systems. Suppose there is a potentially infinite set of paths $x\left(
t\right) _{0}^{T},$ each denoting a particular continuous function $x\left(
t\right) $ for $t\in \left[ 0,T\right] .$ Suppose also that we can evaluate
them, i.e., they give different payoffs. Then, we will learn how to derive
difference or differential equations, that are necessarily satisfied for
optimal paths. If we can solve these equations, we can find the optimal
path. We will use two approaches in this course.

\subsubsection{Dynamic Programming (Bellman).}

Suppose we want to find an optimal investment plan in discrete time and let $%
x_{t}$ denote the stock of capital at time $t.$ Also, let $u_{t}$ denote
investments and assume 
\begin{equation}
x_{t+1}=g(x_{t},u_{t}),
\end{equation}%
which we call the law-of-motion for $x_{t}.$

Each period, the payoff is given by $F\left( t,x_{t},u_{t}\right) $ and the
problem is to solve 
\begin{eqnarray}
&&\max_{x_{0}^{T+1},u_{0}^{T}}\sum_{t=0}^{T}F\left( t,x_{t},u_{t}\right) , \\
\text{ s.t. }x_{t+1} &=&g(x_{t},u_{t})\forall t, \\
x_{0} &=&x.
\end{eqnarray}

Note that this is a dynamic problem; The choice of investment at time $t,$ $%
u_{t}$ may affect payoffs in many future period. First, the payoff at $t,$ $%
F\left( t,x_{t},u_{t}\right) $ is affected directly and also the next
periods payoff since 
\begin{equation*}
F\left( t+1,x_{t+1},u_{t+1}\right) =F\left(
t+1,g(x_{t},u_{t}),u_{t+1}\right) .
\end{equation*}

Furthermore, also payoffs further away can be affected since, for example, $%
x_{t+2}=$ $g(x_{t+1},u_{t+1})=g(g(x_{t},u_{t}),u_{t+1})$, affecting the
payoff in period $t+2.$ The choice of $u_{t}$ must thus take into account
all future payoff relevant effects.

Sometimes the dynamic problem degenerates in the sense that this dynamic
link breaks. To illustrate this, let us simplify and take the example where $%
F\left( t,x_{t},u_{t}\right) =\left( \frac{1}{1+r}\right) ^{t}\left(
f_{t}\left( x_{t}\right) -u_{t}\right) $ and $g\left( x_{t},u_{t}\right)
=\left( 1-\delta \right) x_{t}-u_{t}.$ We can interpret $\delta $ as the
rate of capital depreciation. If $\delta =1$, we have $x_{t+1}=u_{t}$ and by
substituting from the law-of-motion, we can write the problem as%
\begin{eqnarray}
&&\max_{x_{1}^{T}}\sum_{t=0}^{T}\left( \frac{1}{1+r}\right) ^{t}\left(
f_{t}\left( x_{t}\right) -x_{t+1}\right) , \\
x_{0} &=&x.
\end{eqnarray}

The first-order condition for choosing $x_{t}$ for any $t>0$ is $%
f_{t}^{\prime }\left( x_{s}\right) =1+r$, so to know how much to invest in
period $t,$ we only need to know the marginal productivity of capital next
period, i.e., we maximize period-by-period. With $\delta <1,$ we cannot do
this. Then dynamic programing is a handy way to attack the problem. We will
use the Bellman's principle of optimality, saying that there exists a
sequence of functions $V_{t}\left( x_{t}\right) $ such that 
\begin{equation}
V_{t}\left( x_{t}\right) =\max_{u_{t}}F\left( t,x_{t},u_{t}\right)
+V_{t+1}\left( g(x_{t},u_{t})\right)  \label{Bellman}
\end{equation}%
and where%
\begin{eqnarray}
V_{t}\left( x_{t}\right) &\equiv
&\max_{x_{t}^{T+1},u_{t}^{T}}\sum_{s=t}^{T}F\left( s,x_{s},u_{s}\right) \\
\text{ s.t. }x_{t+1} &=&g(x_{t},u_{t}) \\
&&x_{t}\text{ given.}
\end{eqnarray}

$V_{t}\left( x_{t}\right) $ is called the value function and $x_{t}$ a state
variable. We find that through the Bellman principle, we have split up the
problem into a sequence of one period problems. We can then solve the
problem more easily since the first-order condition for maximizing the RHS
of (\ref{Bellman}) 
\begin{equation}
F_{u}\left( t,x_{t},u_{t}\right) +V_{t+1}^{\prime }\left(
g(x_{t},u_{t})\right) g_{u}(x_{t},u_{t})=0.
\end{equation}
implicitly yields a system of difference equations in $x_{t}$ and $u_{t}$
that we may be able to solve.

\subsubsection{Optimal control (Pontryagin)}

The other way of solving dynamic optimization problems that we will use is
called Optimal Control. We will use it for continuous time problems. Suppose
we want to solve 
\begin{eqnarray}
&&\max_{x\left( t\right) _{0}^{T},u\left( t\right)
_{0}^{T}}\int_{0}^{T}F\left( t,x\left( t\right) ,u\left( t\right) \right) dt
\\
s.t.\,\dot{x}\left( t\right) &=&g\left( x\left( t\right) ,u\left( t\right)
\right) \\
x\left( 0\right) &=&x_{0} \\
x\left( T\right) &=&x_{T}.
\end{eqnarray}

Then, Pontryagin's maximum principle says that for each point in time, the
optimal control, call it $u^{\ast }\left( t\right) ,$ satisfies 
\begin{equation}
u^{\ast }\left( t\right) =\arg \max_{u\left( t\right) }F\left( t,x\left(
t\right) ,u\left( t\right) \right) +\lambda \left( t\right) g\left( x\left(
t\right) ,u\left( t\right) \right) ,  \label{Pontryagin}
\end{equation}
where $\lambda \left( t\right) $ can be interpreted as the shadow value of
the state variable (capital). The sum $F\left( t,x\left( t\right) ,u\left(
t\right) \right) +\lambda \left( t\right) g\left( x\left( t\right) ,u\left(
t\right) \right) $, is called the Hamiltonian. Again, we have turned the
dynamic problem into a sequence of static problems. The first order
condition for (\ref{Pontryagin}) will implicitly define a system of
differential equations. Note the similarity between (\ref{Pontryagin}) and (%
\ref{Bellman}).\newpage

\section{Some basics}

\subsection{Taylor series}

We will often need to approximate a continuous and differentiable function
around some point. Specifically, suppose we know $f\left( x_{0}\right) $ and
some of its derivatives and want to approximate the value in some other
point $x\neq x_{0}$. An efficient way of doing such an approximation is to
use Taylors formula. This can be seen as an attempt to fit a polynomial (we
will talk about such below) to a curve.

\begin{result}
The Taylor approximation is%
\begin{eqnarray}
f\left( x\right) &\approx &f\left( x_{0}\right) +\frac{1}{1!}f^{\prime
}\left( x_{0}\right) \left( x-x_{0}\right) +\frac{1}{2!}f^{\prime \prime
}\left( x_{0}\right) \left( x-x_{0}\right) ^{2} \\
&&+\frac{1}{3!}f^{\prime \prime \prime }\left( x_{0}\right) \left(
x-x_{0}\right) ^{3}+...  \notag
\end{eqnarray}
\end{result}

Usually we will just use the first order approximation $f\left( x\right)
\approx f\left( x_{0}\right) +f^{\prime }\left( x_{0}\right) \left(
x-x_{0}\right) $ but sometimes a higher order approximation can be useful.
It can furthermore be shown, that the approximation error of an n'th order
Taylor approximation is given by 
\begin{equation}
\frac{1}{\left( n+1\right) !}\frac{d^{\left( n+1\right) }f\left( c\right) }{%
dx^{\left( n+1\right) }}\left( x-x_{0}\right) ^{n+1}
\end{equation}%
where $c\in \lbrack x_{0},x].$ As we see, when $x-x_{0}$ is small, two
forces imply that the approximation error tends to be small when $n$ is
large. In the denominator, $\left( n+1\right) !$ is large and $\left(
x-x_{0}\right) ^{n+1}$ is small. Although higher-order approximations often
are better, we cannot say that for any $\left( x-x_{0}\right) $ a higher
order approximation is better.

\subsection{Integration\label{integration}}

%TCIMACRO{\TeXButton{TeX field}{\setcounter{equation}{0}}}%
%BeginExpansion
\setcounter{equation}{0}%
%EndExpansion
If $b>a$, the expression 
\begin{equation}
\int\limits_{a}^{b}f\left( t\right) dt,  \label{integral1}
\end{equation}%
can be interpreted as the area under the graph $y=f(x)$ for $x\in \lbrack
a,b]$. How should one compute such an area? The most natural way would be to
divide the interval $[a,b]$ into (many) sub-intervals by choosing numbers $%
a=t_{1}<t_{2}<t_{3}...t_{n}=b$. Suppose we choose uniform intervals, $%
t_{i+1}-t_{i}\equiv \Delta t=\frac{b-a}{n-1}$. Then we can approximate the
area by 
\begin{equation}
\int\limits_{a}^{b}f\left( t\right) dt\approx \sum_{i=1}^{n-1}f\left(
t_{i}\right) \Delta t,  \label{integralapprox}
\end{equation}%
i.e., by summing rectangles of base $\Delta t$ and height $f\left(
t_{i}\right) .$

\FRAME{itbpF}{4.0421in}{1.5255in}{0in}{}{}{slide1.emf}{\special{language
"Scientific Word";type "GRAPHIC";maintain-aspect-ratio TRUE;display
"USEDEF";valid_file "F";width 4.0421in;height 1.5255in;depth
0in;original-width 10.0016in;original-height 7.5135in;cropleft "0";croptop
"1";cropright "1";cropbottom "0.5037";filename
'Figures/Slide1.EMF';file-properties "XNPEU";}}If the function $f\left(
t\right) $ is bounded and differentiable and we let the number of
sub-intervals ($n$) increase and therefore size of each one decrease, the
approximation (\ref{integralapprox}) becomes perfect as $n\rightarrow \infty 
$. To see this, note that we can approximate the error in the approximation
by the triangle $\frac{f\left( t_{i+1}\right) -f\left( t_{i}\right) }{2}%
\Delta t$. Since $f\left( t_{i+1}\right) \approx f\left( t_{i}\right)
+f^{\prime }\left( t_{i}\right) \Delta t$, by a first order Taylor
approximation, each triangle can be approximated by 
\begin{equation}
\frac{f^{\prime }\left( t_{i}\right) \Delta t}{2}\Delta t=\frac{f^{\prime
}\left( t_{i}\right) }{2}\Delta t^{2}.
\end{equation}

Furthermore, the sum of absolute errors, 
\begin{eqnarray}
\left\vert \sum_{i=1}^{n-1}\frac{f^{\prime }\left( t_{i}\right) }{2}\Delta
t^{2}\right\vert &\leq &\max_{i}\left\vert f^{\prime }\left( t_{i}\right)
\right\vert \frac{1}{2}\sum_{i=1}^{n-1}\Delta t^{2}= \\
\max_{i}\left\vert f^{\prime }\left( t_{i}\right) \right\vert \frac{1}{2}%
\left( n-1\right) \left( \frac{b-a}{n-1}\right) ^{2} &=&\max_{i}\left\vert
f^{\prime }\left( t_{i}\right) \right\vert \frac{1}{2}\frac{\left(
b-a\right) ^{2}}{n-1}.
\end{eqnarray}

Clearly, this goes to zero as $n$ goes to infinity. If the function $f\left(
t\right) $ has discontinuities or is non-differentiable somewhere, we can do
the summation for each interval where $f\left( t\right) $ is continuous.

We conclude that we should think of the integral as a sum of rectangles,
each with height given by the function we are integrating and base $dt.$

If we can handle integration over compact intervals, we can also define
integrals over unbounded intervals by taking the limit value as integration
limits ($a$, and/or $b$) approach infinity. Of course, this limit does not
always exist. If it does, we use the notation 
\begin{eqnarray}
\lim_{a\rightarrow -\infty }\int\limits_{a}^{b}f\left( t\right) dt &\equiv
&\int\limits_{-\infty }^{b}f\left( t\right) dt \\
\lim_{b\rightarrow \infty }\int\limits_{a}^{b}f\left( t\right) dt &\equiv
&\int\limits_{a}^{\infty }f\left( t\right) dt
\end{eqnarray}

\subsubsection{Fundamental theorem of calculus}

An integral is a generalization of a sum, and a derivative is a
generalization of a difference. The following theorem links these concepts.
The first part of the\emph{\ Fundamental Theorem} says that if 
\begin{equation}
F\left( b\right) =\int\limits_{a}^{b}f\left( t\right) dt
\end{equation}

then 
\begin{equation}
F^{\prime }\left( b\right) =\frac{\partial }{\partial b}\int\limits_{a}^{b}f%
\left( t\right) dt=f(b).
\end{equation}

(Convince yourself that this is reasonable by making a drawing.) A function $%
F(t)$ that has the property that $F^{\prime }(t)\equiv f(t)$ is called a
primitive for $f$. Clearly, if $F(t)$ is a primitive for $f(t)$, then also $%
F(t)$ plus any constant, is a primitive. Thus, the first part of the
fundamental theorem provides a necessary, but not sufficient condition for
finding the area.

Now let us turn to he second part of the theorem, which provides a way of
calculating the exact value of an integral like in (\ref{integral1}).

\begin{result}
Let $F$ be \emph{any} primitive of $f$, then 
\begin{equation}
\int\limits_{a}^{b}f\left( t\right) dt=F\left( b\right) -F\left( a\right)
\equiv \left[ F\left( t\right) \right] _{a}^{b}.
\end{equation}
\end{result}

Sometimes it is easy to find the primitive, like in the cases $%
f(t)=t^{a},1/(at+b)$ or $e^{at}$, in which case $F(t)$ is $t^{a+1}/(a+1),\ln
(at+b)/a,$ or $e^{at}/a$ respectively. In other cases like $\pi
^{0.5}e^{-x^{2}}$, the primitive cannot be expressed by using standard
algebraic functions. This does not mean that the primitive does not exist.
In the case $\pi ^{0.5}e^{-x^{2}}$, a primitive is the cumulative normal
distribution, which certainly exists.

Since the integral is a (kind of) sum, it is straightforward to understand
Liebniz' rule.

\begin{result}
If $f$ is differentiable, then 
\begin{eqnarray}
\frac{\partial }{\partial x}\int\limits_{a}^{b}f\left( x,t\right) dt
&=&\int\limits_{a}^{b}\frac{\partial }{\partial x}f\left( x,t\right) dt, \\
\frac{\partial }{\partial a}\int\limits_{a}^{b}f\left( x,t\right) dt
&=&-f\left( x,a\right) , \\
\frac{\partial }{\partial b}\int\limits_{a}^{b}f\left( x,t\right) dt
&=&f\left( x,b\right) .
\end{eqnarray}
\end{result}

\subsubsection{Change of variables}

Suppose that $y=g(x)$, then the rules of differentiation gives $dy=g^{\prime
}(x)dx$. Now let us calculate the area under some function $f(y)$ but
integrating over $x.$ In a sense, this is like changing the scale of the
horisontal axis. To do this, we simply substitute $y=g(x),dy=g^{\prime
}(x)dx.$ To get the integration limits, we note that if $x=a$ then $%
y=g\left( a\right) .$

\begin{result}
If, $y=g(x)$ then 
\begin{equation}
\int\limits_{g\left( a\right) }^{g\left( b\right) }f\left( y\right)
dy=\int\limits_{a}^{b}f\left( g\left( x\right) \right) g^{\prime }(x)dx.
\end{equation}
\end{result}

Sometimes a variable substitution makes integration simpler. Take the
following example;%
\begin{equation}
\int\limits_{1}^{2}\left( x^{2}+1\right) ^{10}2xdx.
\end{equation}

Now define $y\equiv x^{2}+1$ $\equiv g\left( x\right) $ implying 
\begin{eqnarray}
dy &=&2xdx,g\left( 1\right) =2,g\left( 2\right) =5 \\
\int\limits_{1}^{2}\left( x^{2}+1\right) ^{10}2xdx
&=&\int\limits_{2}^{5}y^{10}dy
\end{eqnarray}

\subsubsection{Integration by parts}

Another important result, which we will use a lot, is the following rule for
integration by parts. Assume we want to integrate a product of two functions
of $x,$ i.e., $u(x)v(x)$. Then, let $U(x)$ be a primitive of $u$. The rule
for differentiation of products says

\begin{eqnarray}
\frac{d\left( U\left( x\right) v\left( x\right) \right) }{dx} &=&u\left(
x\right) v\left( x\right) +U\left( x\right) v^{\prime }\left( x\right) \\
&\rightarrow &u\left( x\right) v\left( x\right) =\frac{d\left( U\left(
x\right) v\left( x\right) \right) }{dx}-U\left( x\right) v^{\prime }\left(
x\right) .
\end{eqnarray}

Then, we can integrate over $x$ on both sides, giving

\begin{result}
\begin{eqnarray}
\int_{a}^{b}u\left( x\right) v\left( x\right) dx=\int_{a}^{b}\frac{d\left(
U\left( x\right) v\left( x\right) \right) }{dx}dx-\int_{a}^{b}U\left(
x\right) v^{\prime }\left( x\right) dx \\
=\left[ U\left( x\right) v\left( x\right) \right] _{a}^{b}-\int_{a}^{b}U%
\left( x\right) v^{\prime }\left( x\right) dx
\end{eqnarray}
\end{result}

\subsubsection{Double integration}

As we know, the integral is an area under a curve $f(t)$ over an interval $%
[a,b]$ in the $t$ dimension. Similarly, we can compute the volume under a
plane with a height $f(t,s)$ over an area in the $t,s$ dimension. For
example, let $f(t,s)$ be $t\cdot s$ and integrate over a rectangle with
sides $[a_{t},b_{t}]$ and $[a_{s},b_{s}].$

\FRAME{dtbpFU}{3in}{2.0003in}{0pt}{\Qcb{The function, $z=ts$}}{}{Plot}{%
\special{language "Scientific Word";type "MAPLEPLOT";width 3in;height
2.0003in;depth 0pt;display "USEDEF";plot_snapshots TRUE;mustRecompute
FALSE;lastEngine "Maple";xmin "0";xmax "1";ymin "0";ymax "1";xviewmin
"-0.02";xviewmax "1.0204";yviewmin "-0.02";yviewmax "1.0204";zviewmin
"-0.02";zviewmax "1.0204";phi 67;theta -22;plottype 5;plottickdisable
TRUE;axesFont "Times New
Roman,12,0000000000,useDefault,normal";num-x-gridlines 25;num-y-gridlines
25;plotstyle "patchcontour";axesstyle "frame";axestips FALSE;plotshading
"XY";lighting 0;xis \TEXUX{t};yis \TEXUX{s};var1name \TEXUX{$t$};var2name
\TEXUX{$s$};function \TEXUX{$ts$};linestyle 1;pointstyle
"point";linethickness 1;lineAttributes "Solid";var1range "0,1";var2range
"0,1";surfaceColor "[linear:XY:RGB:0x00ff0000:0x000000ff]";surfaceStyle
"Color Patch";num-x-gridlines 25;num-y-gridlines 25;surfaceMesh
"Contour";rangeset"XYZ";valid_file "T";tempfilename
'LDESSQ00.wmf';tempfile-properties "XPR";}}

The volume under the plane $f\left( t,s\right) $ is then given by

\begin{eqnarray}
\int_{a_{s}}^{b_{s}}\int_{a_{t}}^{b_{t}}f\left( t,s\right) dtds
&=&\int_{a_{s}}^{b_{s}}\left( \int_{a_{t}}^{b_{t}}tsdt\right) ds= \\
\int_{a_{s}}^{b_{s}}\left( \left[ \frac{t^{2}}{2}s\right] _{a_{t}}^{b_{t}}%
\right) ds &=&\int_{a_{s}}^{b_{s}}s\left( \frac{b_{t}^{2}-a_{t}^{2}}{2}%
\right) ds= \\
\left( \frac{b_{t}^{2}-a_{t}^{2}}{2}\right) \left[ \frac{s^{2}}{2}\right]
_{a_{s}}^{b_{s}} &=&\left( \frac{b_{t}^{2}-a_{t}^{2}}{2}\right) \left( \frac{%
b_{s}^{2}-a_{s}^{2}}{2}\right)
\end{eqnarray}

Note, that we are first calculating the \emph{area} under $f(t,s)$ over the
interval $[a_{t},b_{t}]$for each $s$. This area is a function of $s$ which
we then integrate over the interval $[a_{s},b_{s}]$ in the $s$-dimension. If
we integrate over other areas than rectangles, the limits of integration are
not independent. For example, we may integrate over a triangle where $%
b_{t}=s $. Simplify and set $a_{t}=a_{s}=0.$ Then, we have

\begin{eqnarray}
\int_{0}^{b_{s}}\int_{0}^{s}f\left( t,s\right) dtds
&=&\int_{0}^{b_{s}}\left( \int_{0}^{s}tsdt\right) ds=\int_{0}^{b_{s}}\left( 
\left[ \frac{t^{2}}{2}s\right] _{0}^{s}\right) ds= \\
\int_{0}^{b_{s}}s\frac{s^{2}}{2}ds &=&\left[ \frac{s^{4}}{8}\right]
_{0}^{b_{s}}=\frac{b_{s}^{4}}{8}.
\end{eqnarray}

\subsection{Complex numbers and trigonometric functions}

The formula for the solution to a quadratic equation $ax^{2}+bx+c=0$, is 
\begin{equation}
x=-\frac{b}{2a}\pm \frac{\sqrt{\left( b^{2}-4ac\right) }}{2a}.
\end{equation}

If $b^{2}-4ac<0$ . The solution is not in the set of real numbers. The
introduction of complex numbers intended to extend the space of solutions to
accommodate such cases and it turns out that for all numbers in the extended
space, we can always find solutions to such equations. We can think of
complex numbers as two-dimension objects $z=(x,y)$. The first number, $x$,
provides the value in the standard real dimension, while the second provides
the value in the other dimension, called \emph{imaginary}. Thus, real
numbers are a special sub-set of complex number such that $y=0.$

\FRAME{itbpF}{3.7732in}{2.5002in}{0in}{}{}{slide3.emf}{\special{language
"Scientific Word";type "GRAPHIC";maintain-aspect-ratio TRUE;display
"USEDEF";valid_file "F";width 3.7732in;height 2.5002in;depth
0in;original-width 10.1918in;original-height 7.5135in;cropleft "0";croptop
"1";cropright "0.5543";cropbottom "0.5037";filename
'Figures/Slide3.EMF';file-properties "XNPEU";}}

The rules for addition and multiplication with complex numbers are the
following

\begin{eqnarray}
z_{1}+z_{2} &=&\left( x_{1}+x_{2},y_{1}+y_{2}\right) \\
z_{1}z_{2} &=&\left( x_{1}x_{2}-y_{1}y_{2},x_{1}y_{2}+x_{2}y_{1}\right)
\end{eqnarray}

Using these rules, we can compute the square of a complex number only
consisting of a unitary imaginary part, i.e., $z=i=\left( 0,1\right) .$

\begin{equation}
\left( 0,1\right) ^{2}=\left( 0-1,0+0\right) =\left( -1,0\right) =-1
\end{equation}

We thus established the important result

\begin{result}
$x=i$ is a solution to the equation $x^{2}=-1$.
\end{result}

Using the rules above, it is also straightforward to show that an
alternative we of writing $z$ is given by the following 
\begin{equation}
z\equiv \left( x,y\right) =\left( x,0\right) +\left( 0,y\right) =x\left(
1,0\right) +y\left( 0,1\right) =x+yi.
\end{equation}

We should also note that

\begin{equation}
\left( x,y\right) \left( x,-y\right) =\left( x^{2}+y^{2},-xy+xy\right)
=\left( x^{2}+y^{2},0\right) \equiv \left| \left( x,y\right) \right| ^{2}
\end{equation}

The numbers $(x,y)$ and $(x,-y)$ are called complex conjugates and the value 
$|x,y|$ is called the modulus of $(x,y)$.

\subsubsection{Polar representation}

Recall that in the right-angled triangle in the figure, we have%
\begin{eqnarray}
\cos \left( \theta \right) &=&\frac{x}{r} \\
\sin \left( \theta \right) &=&\frac{y}{r} \\
r &=&\sqrt{x^{2}+y^{2}}.
\end{eqnarray}

\FRAME{itbpF}{2.6437in}{1.9986in}{0in}{}{}{slide4.emf}{\special{language
"Scientific Word";type "GRAPHIC";maintain-aspect-ratio TRUE;display
"USEDEF";valid_file "F";width 2.6437in;height 1.9986in;depth
0in;original-width 10.0016in;original-height 7.5135in;cropleft "0";croptop
"1";cropright "0.3945";cropbottom "0.6044";filename
'Figures/Slide4.EMF';file-properties "XNPEU";}}

Using this, we can represent the complex number $z$ either by its
coordinates, $(x,y)$ or alternatively as $r(\cos (\theta )+i\sin (\theta ))$%
. The latter form is called the polar representation and $r$ is the \emph{%
modulus} and $\theta $ is the \emph{argument} of $z$. Usually we measure $%
\theta $ in radians.

We will use the following result below.

\begin{result}
Let $z$ be the complex number $(x,y)$ , then $z$ satisfies 
\begin{equation}
z=re^{i\theta }
\end{equation}

where $r$ is the modulus of $z$ and $\theta $ is the argument. In
particular, for $r=1$ we get, 
\begin{eqnarray}
e^{i\theta }=\cos \left( \theta \right) +i\sin \left( \theta \right) \\
e^{-i\theta }=\cos \left( -\theta \right) +i\sin \left( -\theta \right)
=\cos \left( \theta \right) -i\sin \left( \theta \right) \\
e^{i\pi }=\cos \left( \pi \right) +i\sin \left( \pi \right) =-1
\end{eqnarray}
\end{result}

\textbf{Optional proof:}

The Taylor formula around zero for any function $f\left( x\right) $ is 
\begin{equation}
f\left( x\right) =f\left( 0\right) +\frac{f^{\prime }\left( 0\right) }{1!}x+%
\frac{f^{\prime \prime }\left( 0\right) }{2!}x^{2}+\frac{f^{\prime \prime
\prime }\left( 0\right) }{3!}x^{3}...
\end{equation}

Using this for $f(x)=e^{x},\cos (x)$ and $\sin (x)$, and the rules 
\begin{eqnarray}
\frac{\partial }{\partial x}\cos \left( x\right) &=&-\sin \left( x\right) \\
\frac{\partial }{\partial x}\sin \left( x\right) &=&\cos \left( x\right)
\end{eqnarray}%
we have respectively 
\begin{eqnarray}
e^{x} &=&1+x+\frac{x^{2}}{2!}+\frac{x^{3}}{3!}+\frac{x^{4}}{4!}... \\
\cos \left( x\right) &=&\cos \left( 0\right) -\sin (0)x-\frac{\cos \left(
0\right) }{2!}x^{2}+\frac{\sin \left( 0\right) }{3!}x^{3} \\
&&+\frac{\cos \left( 0\right) }{4!}x^{4}-\frac{\sin \left( 0\right) }{5!}%
x^{5}-\frac{\cos \left( 0\right) }{6!}x^{6}... \\
&=&1-\frac{x^{2}}{2}+\frac{x^{4}}{4!}-\frac{x^{6}}{6!}... \\
\sin \left( x\right) &=&\sin \left( 0\right) +\cos (0)x-\frac{\sin \left(
0\right) }{2!}x^{2}-\frac{\cos \left( 0\right) }{3!}x^{3} \\
&&+\frac{\sin \left( 0\right) }{4!}x^{4}+\frac{\cos \left( 0\right) }{5!}%
x^{5}-\frac{\sin \left( 0\right) }{6!}x^{6}-\frac{\cos \left( 0\right) }{7!}%
x^{7}... \\
&=&x-\frac{1}{3!}x^{3}+\frac{1}{5!}x^{5}-\frac{1}{7!}x^{7}
\end{eqnarray}

Thus, using the Taylor formula around zero to evaluate $f(i\theta
)=e^{i\theta }$, we have

\begin{eqnarray}
e^{i\theta } &=&1+i\theta +\frac{\left( i\theta \right) ^{2}}{2!}+\frac{%
\left( i\theta \right) ^{3}}{3!}+\frac{\left( i\theta \right) ^{4}}{4!}+%
\frac{\left( i\theta \right) ^{5}}{5!}+\frac{\left( i\theta \right) ^{6}}{6!}%
... \\
&=&1+i\theta -\frac{1}{2}\theta ^{2}-\frac{1}{3!}i\theta ^{3}+\frac{1}{4!}%
\theta ^{4}+\frac{1}{5!}i\theta ^{5}-\frac{1}{6!}\theta ^{6}... \\
&=&1-\frac{1}{2}\theta ^{2}+\frac{1}{4!}\theta ^{4}-\frac{1}{6!}\theta
^{6}... \\
&&+i\theta -\frac{1}{3!}i\theta ^{3}+\frac{1}{5!}i\theta ^{5}... \\
&=&\cos \left( \theta \right) +i\sin \left( \theta \right)
\end{eqnarray}

\subsubsection{Polynomials}

A polynomial $P(z)$ of order $n$ is defined as weighted sum of $z^{s}$ for $%
s\in \left\{ 0,...,n\right\} $, i.e., 
\begin{equation}
P(z)\equiv a_{n}z^{n}+a_{n-1}z^{n-1}+,...,+a_{1}z^{1}+a_{0},
\end{equation}
for some sequence of constants $\left\{ a_{s}\right\} _{0}^{n}.$ The
following will be important for our analysis of difference and differential
equations.

\begin{result}
A polynomial $P(z)$ of order $n$ has exactly $n$, not necessarily distinct,
roots. I.e., it can be expressed as
\end{result}

\begin{equation}
P(z)=a_{n}\left( z-r_{1}\right) \left( z-r_{2}\right) ...\left(
z-r_{n}\right)  \label{Poly1}
\end{equation}

From (\ref{Poly1}), we see that each root $r_{i}$ is a solution to the
equation $P(z)=0$. Note that the roots may be repeated, i.e.,$r_{i}=r_{j}$
and that, of course, some root may be complex. It turns out, also, that
complex roots always come in pairs, the complex conjugates.

\newpage

\section{Differential equations}

\subsection{Linear Differential Equations of First Order}

%TCIMACRO{\TeXButton{TeX field}{\setcounter{equation}{0}}}%
%BeginExpansion
\setcounter{equation}{0}%
%EndExpansion
A first order differential equation is an equation of the form 
\begin{equation}
\frac{dx(t)}{dt}\equiv \dot{x}\left( t\right) =f(x,t).  \label{diffeq1}
\end{equation}%
As noted above, there will in general be a whole class of functions $x(t;c)$
(parameterized by $c$) that satisfies the differential equation (\ref%
{diffeq1}). We need more information, like an initial condition for $x\left(
t_{0}\right) ,$ to pin down the solution exactly.

\begin{result}
Given that $f$ is continuous and has continuous first derivatives, there is
going to be a one function $x(t)$ that satisfies (\ref{diffeq1}) and an
initial condition.
\end{result}

\subsubsection{The simplest case}

If $f$ is independent of $x$, the solution is trivial. Then since%
\begin{equation}
\dot{x}\left( t\right) =f\left( t\right) ,  \label{simplestcase}
\end{equation}%
the class of functions satisfying this is a primitive function of $f$ plus
any constant, i.e., for any $t_{0}$ and $\tilde{c},$ 
\begin{equation}
x\left( t\right) =\int_{t_{0}}^{t}f\left( s\right) ds+\tilde{c}=F\left(
t\right) -F\left( t_{0}\right) +\tilde{c},  \label{solutiondiff1b}
\end{equation}%
where $F\left( t\right) $ is any primitive of $f$, satisfies (\ref%
{simplestcase}). Note that $F(t_{0})$ is a constant. We can thus merge the
two constants, defining $c\equiv $ $\tilde{c}-F\left( 0\right) $ and write 
\begin{equation}
x\left( t\right) =F\left( t\right) +c,  \label{simplestcasesol}
\end{equation}%
which shows that there is only one degree of freedom in the constants $\hat{c%
}$ and $t_{0}$. Choosing another $t_{0}$ simply means that the constant $%
\hat{c}$ has to be chosen in another way. Certainly, for all $c$, (\ref%
{simplestcasesol}) satisfies (\ref{diffeq1}). For example, if 
\begin{eqnarray}
f\left( t\right) &=&e^{at}, \\
F(t) &=&\frac{e^{at}}{a}
\end{eqnarray}%
is a primitive for $f\left( t\right) .$ The arbitrary constant is pinned
down with some other piece of information. So, if we want to find $x(t)$ and
we know the value of $x(0)$, we get 
\begin{eqnarray}
x\left( t\right) &=&F\left( t\right) +c \\
x\left( 0\right) &=&F\left( 0\right) +c \\
&\rightarrow &c=x\left( 0\right) -F\left( 0\right) \\
x\left( t\right) &=&F\left( t\right) +x\left( 0\right) -F\left( 0\right) .
\end{eqnarray}

\subsubsection{A note on notation}

Above, we have used the proper notation of an integral, where both the lower
and upper limits and the dummy variable to integrate over have separate
names and are all written out. Often, a more sloppy notation is used, for an
arbitrary $t_{0}$ we write%
\begin{equation}
\int_{t_{0}}^{t}f\left( s\right) ds=\int f\left( t\right) dt=F\left( t\right)
\end{equation}%
where it is understood that $\int f\left( t\right) dt$ is a (any) primitive
of $f\left( t\right) .$ This notation, called an indefinite integral, saves
on variables, but can be confusing since $t$ is used as both the variable to
integrate over and as upper integration limit. Nevertheless, I will follow
ordinary practice and use it below. Using this notation, we rewrite (\ref%
{solutiondiff1b})

\begin{equation}
x\left( t\right) =\int f\left( t\right) dt+c.
\end{equation}

\subsubsection{A little bit more complicated}

Very often, a solution to a more complicated differential equation is
derived by transforming the original equation into something that has the
form of (\ref{simplestcase}). Linear first order differential with constant
coefficients equations can be solved directly using such a transformation.
Consider 
\begin{equation}
\dot{y}\left( t\right) +ry\left( t\right) =q.  \label{arbitrage1}
\end{equation}

In this case, we multiply both sides by $e^{rt}$ (often called the
integrating factor). After doing that, note that the LHS becomes 
\begin{equation}
e^{rt}\left( \dot{y}\left( t\right) +ry\left( t\right) \right) =\frac{%
d\left( e^{rt}y\left( t\right) \right) }{dt}.
\end{equation}

Thus, thinking of $e^{rt}y\left( t\right) $ as simply a function of $t$, as $%
x(t)$ in (\ref{simplestcase}), we get a LHS that is the time derivative of a
known function of $t$ and the RHS is also a function \emph{only} of $t$.
Let's set $t_{0}=0,$ then the solution is found as in (\ref{solutiondiff1b}%
). 
\begin{eqnarray}
\frac{d\left( e^{rt}y\left( t\right) \right) }{dt} &=&e^{rt}q \\
e^{rt}y\left( t\right) &=&\int_{0}^{t}e^{rs}qds+c=\left[ \frac{e^{rs}q}{r}%
\right] _{0}^{t}+c \\
&=&\frac{q}{r}\left( e^{rt}-1\right) +c \\
&\rightarrow &y\left( t\right) =e^{-rt}\left( \frac{q}{r}\left(
e^{rt}-1\right) +c\right) =\frac{q}{r}+\left( c-\frac{q}{r}\right) e^{-rt}.
\end{eqnarray}

If we know $y(t)$ at some point in time, e.g. at$\ t=0.$%
\begin{eqnarray}
y\left( 0\right) &=&\frac{q}{r}+\left( c-\frac{q}{r}\right) =c \\
&\rightarrow &y\left( t\right) =\frac{q}{r}+\left( y\left( 0\right) -\frac{q%
}{r}\right) e^{-rt}
\end{eqnarray}

Let's compare this to the case when we use the indefinite integral. Then, we
have 
\begin{eqnarray}
\frac{d\left( e^{rt}y\left( t\right) \right) }{dt} &=&e^{rt}q \\
e^{rt}y\left( t\right) &=&\int e^{rt}qdt+c=\frac{e^{rt}q}{r}+c \\
&=&\frac{q}{r}e^{rt}+c \\
&\rightarrow &y\left( t\right) =\frac{q}{r}+ce^{-rt}.
\end{eqnarray}%
\begin{eqnarray}
y\left( 0\right) &=&\frac{q}{r}+c\rightarrow c=\left( y\left( 0\right) -%
\frac{q}{r}\right) \\
&\rightarrow &y\left( t\right) =\frac{q}{r}+\left( y\left( 0\right) -\frac{q%
}{r}\right) e^{-rt}
\end{eqnarray}

As you see, although the integration constants are different, the result is
the same. The differerence is that the constant submerges the primitive
evaluated at the lower integration limit.

\subsubsection{Variable exogenous parts}

Let us generalize (\ref{arbitrage1}) by assuming that $q=q\left( t\right) $,
a function of time. Consider an arbitrage condition where $y(t)$ is an asset
price, with dividends $q\left( t\right) $ per unit of time$.$ Not to loose
track of things, we use the notation with integration limits, rather than
indefinite integrals.

Holding the asset should yield no excess return, i.e., 
\begin{eqnarray}
\dot{y}\left( t\right) +q\left( t\right) &=&ry\left( t\right) \\
\dot{y}\left( t\right) -ry\left( t\right) &=&-q\left( t\right) \\
&\Rightarrow &e^{-rt}y\left( t\right) =-\int_{t_{0}}^{t}e^{-rs}q\left(
s\right) ds+c
\end{eqnarray}

Often, economic theory let us use a so called \textquotedblright
No-Ponzi\textquotedblright\ condition (typically requiring $r>0)$ . Here
this gives an end condition, since it requires that 
\begin{equation}
\lim_{t\rightarrow \infty }e^{-rt}y\left( t\right) =0.
\end{equation}%
Using this, we get 
\begin{equation}
c=\lim_{t\rightarrow \infty }\int_{t_{0}}^{t}e^{-rs}q\left( s\right)
ds\equiv \int_{t_{0}}^{\infty }e^{-rs}q\left( s\right) ds,
\end{equation}%
giving 
\begin{eqnarray}
e^{-rt}y\left( t\right) &=&\int_{t_{0}}^{\infty }e^{-rs}q\left( s\right)
d-\int_{t_{0}}^{t}e^{-rs}q\left( s\right) ds \\
&=&\int_{t}^{\infty }e^{-rs}q\left( s\right) ds \\
y\left( t\right) &=&e^{rt}\int_{t}^{\infty }e^{-rs}q\left( s\right)
ds=\int_{t}^{\infty }e^{-r\left( s-t\right) }q\left( s\right) ds,
\end{eqnarray}%
i.e., that no arbitrage and \textquotedblright No-Ponzi\textquotedblright\
implies that the price must equal the discounted sum of future dividends.

\subsubsection{Variable coefficients}

If also the coefficient on the endogenous variable ($y\left( t\right) $) is
varying over time, we need a more general integrating factor to make the LHS
od the equation being the time differential of a known function. For
example, consider the no-arbitrage equation for an asset price $y\left(
t\right) $ with dividends $q\left( t\right) $ and interest rate $r\left(
t\right) $%
\begin{equation}
\dot{y}\left( t\right) +q\left( t\right) =r\left( t\right) y\left( t\right)
\Rightarrow \dot{y}\left( t\right) -r\left( t\right) y\left( t\right)
=-q\left( t\right) .
\end{equation}

Here the integrating factor is 
\begin{equation}
e^{-\int_{t_{0}}^{t}r(s)ds}.
\end{equation}%
with 
\begin{equation}
\frac{de^{-\int_{t_{0}}^{t}r(s)ds}}{dt}=-r\left( t\right)
e^{-\int_{t_{0}}^{t}r(s)ds},
\end{equation}%
alternatively expressed as%
\begin{equation*}
\frac{de^{^{-\int r(t)dt}}}{dt}=-r\left( t\right) e^{-\int r(t)dt}.
\end{equation*}

Approximating the integral as a sum of rectangles with base $\Delta t$ as in
section \ref{integration}, (which is exact if $r(t)$ were piecewise
constant), and defining $r\left( t_{0}+s\Delta t\right) \equiv r_{t_{s}}$,
for the integers $s\in \left\{ 0,S\right\} ,$ $S=(t-t_{0})/\Delta t$, the
integrating factor, can be written 
\begin{equation}
e^{-\int_{t_{0}}^{t}r(s)ds}\approx e^{-r_{0}\Delta t}e^{-r_{1}\Delta
t}...e^{-r_{S}\Delta t}\approx \left( \frac{1}{1+r_{1}\Delta t}\right)
\left( \frac{1}{1+r_{1}\Delta t}\right) ...\left( \frac{1}{1+r_{S}\Delta t}%
\right) ,
\end{equation}%
i.e., it is product of all short run discount factors between $t_{0}$ and $%
t. $ To save on notation, this product is denoted as 
\begin{equation}
e^{-\int_{t_{0}}^{t}r(s)ds}\equiv R\left( t;t_{0}\right) ,
\end{equation}%
or if the starting point is suppressed as $R\left( t\right) .$ This has a
straightforward interpretation. Suppose the variable discount rate is given
by $r\left( s\right) $, then, the the discounted value of a payment $y\left(
t\right) $ at $t$ seen from $t_{0}$ is $R\left( t;t_{0}\right) y\left(
t\right) .$Using, the integrating factor, we get 
\begin{eqnarray}
R\left( t\right) \left( \dot{y}\left( t\right) -r\left( t\right) y\left(
t\right) \right) &=&-R\left( t\right) q\left( t\right) \\
\frac{dR\left( t\right) y\left( t\right) }{dt} &=&-R\left( t\right) q\left(
t\right) ,\rightarrow \\
R\left( t\right) y\left( t\right) &=&-\int_{t_{0}}^{t}R\left( s\right)
q\left( s\right) ds+c
\end{eqnarray}

Suppose again, $\lim {}_{t\rightarrow \infty }R\left( t\right) y\left(
t\right) =0,$ implying $c=$ $\int_{t_{0}}^{\infty }R\left( s\right) q\left(
s\right) ds$. Then, noting that $R\left( t\right) ^{-1}R\left( s\right)
=e^{\int_{t_{0}}^{t}r(v)dv-\int_{t_{0}}^{s}r(v)dv}=e^{-\int_{t}^{s}r(v)dv}$ $%
=R\left( s;t\right) $ we have 
\begin{eqnarray}
R\left( t\right) y\left( t\right) &=&-\int_{t_{0}}^{t}R\left( s\right)
q\left( s\right) ds+\int_{t_{0}}^{\infty }R\left( s\right) q\left( s\right)
ds=\int_{t}^{\infty }R\left( s\right) q\left( s\right) ds \\
y\left( t\right) &=&\int_{t}^{\infty }R\left( t\right) ^{-1}R\left( s\right)
q\left( s\right) ds=\int_{t}^{\infty }e^{-\int_{t}^{s}r(v)dv}q\left(
s\right) ds,
\end{eqnarray}%
i.e., that the asset price equals the discounted value of future dividends.

As another example, consider money on a bank account with variable interest
rate and deposits $q\left( t\right) ,$then 
\begin{eqnarray}
\dot{y}\left( t\right) &=&r\left( t\right) y\left( t\right) +q\left( t\right)
\\
\dot{y}\left( t\right) -r\left( t\right) y\left( t\right) &=&q\left( t\right)
\\
\frac{dR\left( t,t_{0}\right) y\left( t\right) }{dt} &=&R\left(
t,t_{0}\right) q\left( t\right) \\
&\rightarrow &R\left( t,t_{0}\right) y\left( t\right)
=\int_{t_{0}}^{t}R\left( s,t_{0}\right) q\left( s\right) ds+c \\
y\left( t\right) &=&\int_{t_{0}}^{t}R\left( t,t_{0}\right) ^{-1}R\left(
s,t_{0}\right) q\left( s\right) +R\left( t,t_{0}\right) ^{-1}c
\end{eqnarray}

Since $s\leq t$ here, the term $R\left( t,t_{0}\right) ^{-1}R\left(
s,t_{0}\right) $ is more conveniently\footnote{%
But remember $-\int_{t}^{s}r(v)dv=\int_{s}^{t}r(v)dv.$} written $%
e^{\int_{s}^{t}r(v)dv}.$ Using also an initial condition, $y\left(
t_{0}\right) ,$we have 
\begin{equation}
y\left( t\right) =\int_{t_{0}}^{t}e^{\int_{s}^{t}r(v)dv}q\left( s\right)
ds+e^{\int_{t_{0}}^{t}r(v)dv}y\left( t_{0}\right) ,
\end{equation}%
where the first term is the period $t$ value of all deposits up until $t$
and the second term is the period $t$ value of the initial amount on the
bank.

\subsubsection{Separating variables}

Sometimes, we can write a differential equation such that the LHS only
contains functions of $x$ and $\dot{x}$ and the RHS only a function of $t$.
For example 
\begin{eqnarray}
\dot{x}\left( t\right) &=&\frac{h\left( t\right) }{g\left( x\right) } \\
\dot{x}\left( t\right) g\left( x\right) &=&h\left( t\right) .
\end{eqnarray}

In this case, we can use the following trick. Let $G\left( x\right) $ be any
primitive of $g\left( x\right) $, i.e., $\frac{dG\left( x\right) }{dx}%
=g\left( x\right) .$ Then,

\begin{eqnarray}
\dot{x}\left( t\right) g\left( x\left( t\right) \right) &=&\frac{dG(x\left(
t\right) )}{dt}=h\left( t\right) \\
G(x) &=&\int h\left( s\right) ds+c
\end{eqnarray}

We can then recover $x$ by inverting $G$. An example;

\begin{eqnarray}
\dot{x}\left( t\right) &=&\left( x\left( t\right) t\right) ^{2} \\
\dot{x}\left( t\right) x\left( t\right) ^{-2} &=&t^{2}.
\end{eqnarray}%
Now let $g\left( x\right) =x^{-2}$ implying $G\left( x\right) =-x^{-1}$.
Then since $\frac{dG\left( x\left( t\right) \right) }{dt}=\dot{x}\left(
t\right) g\left( x\left( t\right) \right) =\dot{x}\left( t\right) x\left(
t\right) ^{-2}$, we have 
\begin{eqnarray}
\frac{dG\left( x\left( t\right) \right) }{dt} &=&t^{2} \\
G\left( x\left( t\right) \right) &=&-x\left( t\right) ^{-1}=\int t^{2}dt+c=%
\frac{t^{3}}{3}+c.
\end{eqnarray}%
Using the fact that $z=G\left( x\right) =-x^{-1}\Rightarrow x=-z^{-1}$, so
the inverse function is given by $G\left( z\right) =-z^{-1}$, we get 
\begin{equation*}
x\left( t\right) =G^{-1}\left( \frac{t^{3}}{3}+c\right) =-\left( \frac{t^{3}%
}{3}+c\right) ^{-1}.
\end{equation*}

\subsection{Linear differential equations of higher order}

\subsubsection{Linear second order differential equations}

A linear second order differential equation has the form 
\begin{equation}
\ddot{y}\left( t\right) +p\left( t\right) \dot{y}\left( t\right) +q\left(
t\right) y\left( t\right) =R\left( t\right)
\end{equation}

This cannot be solved directly by transformations in the simple way we did
with first order equations. Instead we use a more general method (that would
have worked above also). First some definitions

\begin{definition}
The homogeneous part of differential equation is achieved by setting all
exogenous variables, including constants, to zero.
\end{definition}

\begin{definition}
Two functions $y_{1}(t)$and $y_{2}(t)$ are linearly independent in a region $%
\Omega $ iff there is no $c_{1},c_{2}\neq \{0,0\}$ s.t. 
\begin{equation}
c_{1}y_{1}(t)=c_{2}y_{2}(t)\text{ }\forall t\in \Omega .
\end{equation}
\end{definition}

\begin{result}
$\label{solutiondiff1}$The general solution (the complete set of solutions)
to a differential equation is the general solution to the homogeneous part
plus any particular solution to the complete equation.

The general solution of the homogeneous part of a second order linear
differential equation can be expressed as $c_{1}y_{1}(t)+c_{2}y_{2}(t)$,
where $y_{1}(t)$ and $y_{2}(t)$ are two linearly independent particular
solutions to the homogeneous equations and $c_{1}$ and $c_{2}$ are some
arbitrary integration constants.
\end{result}

\subsubsection{Homogeneous Equations with Constant Coefficients}

Consider the homogeneous part of a differential equation, given by 
\begin{equation}
\ddot{y}\left( t\right) +p\dot{y}\left( t\right) +qy\left( t\right) =0.
\end{equation}
To solve this equation, we first define the \emph{characteristic equation, }%
for this equation; 
\begin{equation}
r^{2}+pr+q=0.
\end{equation}

Since this is a second-order polynomial, it has two roots 
\begin{equation}
r_{1,2}=-\frac{1}{2}p\pm \frac{1}{2}\sqrt{\left( p^{2}-4q\right) }.
\end{equation}

Now, it is straightforward to see that $e^{r_{1}t}$ and $e^{r_{2}t}$ both
are solutions to the homogeneous equation, by noting that for $i\in \left\{
1,2\right\} $%
\begin{eqnarray}
\frac{de^{r_{i}t}}{dt} &=&r_{i}e^{r_{i}t} \\
\frac{d^{2}e^{r_{i}t}}{dt^{2}} &=&r_{i}^{2}e^{r_{i}t}\rightarrow \\
\frac{d^{2}e^{r_{i}t}}{dt^{2}}+p\frac{de^{r_{i}t}}{dt}+qe^{r_{i}t} &=&\left(
r_{i}^{2}+pr_{i}+q\right) e^{r_{i}t}=0.
\end{eqnarray}

So then, result \ref{solutiondiff1} tells us that the general solution to
the homogenous equation is 
\begin{equation}
y_{h}\left( t\right) =c_{1}e^{r_{1}t}+c_{2}e^{r_{2}t},
\end{equation}%
\emph{provided }$e^{r_{1}t}$ and $e^{r_{2}t}$ are linearly independent. It
is easy to verify that this is the case, if and only if $r_{1}\neq r_{2}.$

\subsubsection{Complex roots}

Complex roots pose no particular difficulty, we simply have to recall that
for any real number $b$%
\begin{eqnarray}
e^{bi} &=&\cos \left( b\right) +i\sin \left( b\right) , \\
e^{-bi} &=&\cos \left( b\right) -i\sin \left( b\right) ,
\end{eqnarray}
yielding for the complex roots $r_{1,2}=a\pm bi,$%
\begin{gather}
c_{1}e^{\left( a+bi\right) t}+c_{2}e^{\left( a-bi\right) t}= \\
c_{1}e^{at}\left( \cos \left( bt\right) +i\sin \left( bt\right) \right)
+c_{2}e^{at}\left( \cos \left( bt\right) -i\sin \left( bt\right) \right) \\
=e^{at}\left( \left( c_{1}+c_{2}\right) \cos \left( bt\right) +\left(
c_{1}-c_{2}\right) i\sin \left( bt\right) \right)
\end{gather}

Note here, that $c_{1}$ and $c_{2}$ are any constants in the space of
complex numbers are. Defining 
\begin{eqnarray}
c_{1}+c_{2} &\equiv &\bar{c}_{1}, \\
\left( c_{1}-c_{2}\right) i &\equiv &\bar{c}_{2},
\end{eqnarray}

It turns out that for any $\bar{c}_{1}$ and $\bar{c}_{2}$ on the real line,
we can find $c_{1}$ and $c_{2}$ satisfying this definition. Since we, at
least in economics, are only interested in solutions in the real space, we
can use the restricted set of constants satisfying $\bar{c}_{1}$ and $\bar{c}%
_{2}$ being on the real line. We can then write the general solution \emph{%
in the real space} as 
\begin{equation}
y_{h}\left( t\right) =e^{at}\left( \bar{c}_{1}\cos \left( bt\right) +\bar{c}%
_{2}\sin \left( bt\right) \right) .
\end{equation}

\subsubsection{Repeated roots}

The general solution to the homogenous equation in the case when the roots
are repeated, i.e., $r_{1}=r_{2}\equiv r$ is 
\begin{equation}
y_{h}\left( t\right) =c_{1}e^{rt}+c_{2}te^{rt}.
\end{equation}

Convince yourself that they are linearly independent and check that they are
both solutions if the roots are repeated but not otherwise!

\subsubsection{Non-Homogeneous Equations with Constant Coefficients}

Relying on result \ref{solutiondiff1}, the only added problem when we have a
non-homogeneous equation is that we have to find \emph{one} particular
solution to the complete equation. Consider 
\begin{equation}
\ddot{y}\left( t\right) +p\dot{y}\left( t\right) +qy\left( t\right) =R\left(
t\right) .
\end{equation}
Typically we guess a form of this solution and then use the \emph{method of
undetermined coefficients}. Often a good guess is a solution of a form
similar to $R(t)$, e.g., if it is a polynomial of degree $n$ we guess on a
general $n^{\prime }th$ degree polynomial with unknown coefficients. The
simplest example is if $R\left( t\right) $ equals a constant $R$, we then
guess on a constant, $y_{p}\left( t\right) =y^{ss}$, i.e., a \emph{steady
state}, in which case $\ddot{y}\left( t\right) $ and $\dot{y}\left( t\right) 
$ both are zero. To satisfy the differential equation, 
\begin{eqnarray}
qy^{ss} &=&R \\
y^{ss} &=&\frac{R}{q}.
\end{eqnarray}

Another example is 
\begin{equation}
\ddot{y}\left( t\right) -2\dot{y}\left( t\right) +y\left( t\right) =3t^{2}+t,
\end{equation}

in which case we guess 
\begin{equation}
y_{p}\left( t\right) =At^{2}+Bt+C,
\end{equation}

for some, yet \emph{undetermined coefficients }$A,B$ and $C$\emph{. }We then
solve for these constants by substituting into the differential equation 
\begin{equation}
2A-2\left( 2At+B\right) +At^{2}+Bt+C=3t^{2}+t.
\end{equation}

For this to hold for each $t$, we need 
\begin{eqnarray}
A &=&3 \\
-4A+B &=&1 \\
2A-2B+C &=&0
\end{eqnarray}
yielding $A=3,B=13,$ and $C=20.$ So a particular solution is 
\begin{equation}
y_{p}\left( t\right) =3t^{2}+13t+20.
\end{equation}

The characteristic equation is 
\begin{equation}
r^{2}-2r+1=\left( r-1\right) ^{2}\rightarrow r_{1,2}=1.
\end{equation}
So the general solution is 
\begin{equation}
y\left( t\right) =c_{1}e^{t}+c_{2}te^{t}+3t^{2}+13t+20.
\end{equation}

\subsubsection{Linear nth Order Differential Equations}

Consider the nth order differential equation 
\begin{equation}
y^{n}\left( t\right) +P_{1}y^{n-1}\left( t\right) +...+P_{n}y\left( t\right)
=R\left( t\right) ,  \label{norderdiff}
\end{equation}
where 
\begin{equation}
y^{n}\left( t\right) \equiv \frac{d^{n}y\left( t\right) }{dt^{n}}.
\end{equation}

The solution technique is here exactly analogous to the second order case.
First, some definitions

\begin{definition}
For any differentiable function $y\left( t\right) $, the differential
operator $D$ is defined as 
\begin{equation}
Dy\left( t\right) \equiv \frac{dy\left( t\right) }{dt},
\end{equation}
satisfying $D^{n}y\left( t\right) \equiv \frac{d^{n}y\left( t\right) }{dt^{n}%
}.$
\end{definition}

Using this definition, and defining the following polynomial, 
\begin{equation}
P\left( r\right) \equiv r^{n}+P_{1}r^{n-1}+...+P_{n},
\end{equation}%
we have 
\begin{eqnarray}
P\left( D\right) y\left( t\right) &=&\left(
D^{n}+P_{1}D^{n-1}+...+P_{n}\right) y\left( t\right) \\
&=&D^{n}y\left( t\right) +P_{1}D^{n-1}y\left( t\right) +...+P_{n}y\left(
t\right)
\end{eqnarray}%
so we can write (\ref{norderdiff}) in the condensed form, 
\begin{equation}
P\left( D\right) y\left( t\right) =R\left( t\right)
\end{equation}%
and the characteristic equation is 
\begin{equation}
P\left( r\right) =0,
\end{equation}%
with roots $r_{1,...,n}$

The general solution to the homogenous part is then the sum of the $n$
solutions corresponding to each of $n$ roots. The only thing to note is that
if one root, $r$, is repeated $k\geq 2$ times, the solutions corresponding
to this root are given by 
\begin{equation}
c_{1}e^{rt}+c_{2}te^{rt}+...+c_{k}t^{k-1}e^{rt}.
\end{equation}

Repeated complex roots are handled the same way. Say the root $a\pm bi$ is
repeated in $k$ pairs. Their contribution to the general solution is given
by 
\begin{eqnarray}
&&e^{at}\left( c_{1}\cos \left( bt\right) +c_{2}\sin \left( bt\right)
+t\left( c_{3}\cos \left( bt\right) +c_{4}\sin \left( bt\right) \right)
\right) +... \\
&&+e^{at}t^{k-1}\left( c_{2k-1}\cos \left( bt\right) +c_{2k}\sin \left(
bt\right) \right) .
\end{eqnarray}

A particular solution to the complete equations can often be solved by the
method of undetermined coefficients.

\subsection{Stability}

From the solutions to the differential equations we have seen we find that
the terms corresponding to roots that have positive real parts (unstable, or
explosive roots) tend to explode as $t$ goes to infinity. This means that
also the solution explodes unless the corresponding integration constants
are zero. Terms with roots that have negative real parts (stable roots), on
the other hand, always converge to zero.

If all roots are \emph{strictly} negative for a differential equation, $%
P\left( D\right) y\left( t\right) =r$, the system converges to a unique
point as time goes to infinity, wherever it starts. This point is often
called a \emph{sink} and the system is defined as globally stable.

If all roots are strictly positive the system is globally unstable, but
there is still a steady-state, this is sometimes called the \emph{origin}.
The reason for this is that if the system is unstable in all dimensions when
time goes forward, it will be \emph{stable} in all dimensions if time goes
in reverse. Starting from \emph{any} point and going \emph{backward, }one
reaches, in the limit the steady state, i.e., it is the origin of all paths.

If a system has both stable and unstable roots, it is called saddle-path
stable. Then, in some sub-dimensions it is stable.

\subsubsection{Non-linear equations and local stability}

Look at the nonlinear differential equation 
\begin{equation}
\dot{x}\left( t\right) =x\left( t\right) ^{2}-4.
\end{equation}

Although we have not learned how to solve such an equation, we can say
something about it. Let us plot the relation $x\left( t\right) \rightarrow 
\dot{x}\left( t\right) ;$

\FRAME{itbpF}{3.6478in}{3.0026in}{0in}{}{}{slide5.emf}{\special{language
"Scientific Word";type "GRAPHIC";maintain-aspect-ratio TRUE;display
"USEDEF";valid_file "F";width 3.6478in;height 3.0026in;depth
0in;original-width 10.0016in;original-height 7.5135in;cropleft "0";croptop
"1";cropright "0.5458";cropbottom "0.4029";filename
'Figures/Slide5.EMF';file-properties "XNPEU";}}

We see that $x(t)=2$ and $x(t)=-2$ are stationary points. We also see that $%
x(t)=-2$ is locally stable. In the region $[-\infty ,2)$ $x$ converge to $%
x=-2$. 2 is an unstable stationary point and in the region $(2,\infty ]$, $x$
explodes over time. In a plot $x\left( t\right) \rightarrow \dot{x}\left(
t\right) ,$ local stability is equivalent to a negative slope at the
stationary point.

\subsection{Systems of Linear First Order Differential Equations}

Consider the following system of two first order differential equations 
\begin{eqnarray}
\left[ 
\begin{array}{c}
\dot{y}_{1}\left( t\right) \\ 
\dot{y}_{2}\left( t\right)%
\end{array}%
\right] &=&\left[ 
\begin{array}{cc}
a_{11} & a_{12} \\ 
a_{21} & a_{22}%
\end{array}%
\right] \left[ 
\begin{array}{c}
y_{1}\left( t\right) \\ 
y_{2}\left( t\right)%
\end{array}%
\right] +\left[ 
\begin{array}{c}
p_{1}\left( t\right) \\ 
p_{2}\left( t\right)%
\end{array}%
\right]  \label{diffsystem1} \\
\mathbf{\dot{y}}\left( t\right) &=&\mathbf{Ay}\left( t\right) +\mathbf{p}%
\left( t\right)
\end{eqnarray}

As in the one equations case we start by finding the general solutions to
the homogeneous part. This plus some particular solution is the general
solution to the complete system.

If the off diagonal terms are zero the solution to the homogeneous part is
trivial, since there is no interdependency between the equations 
\begin{equation}
\left[ 
\begin{array}{c}
\dot{y}_{1}\left( t\right) \\ 
\dot{y}_{2}\left( t\right)%
\end{array}
\right] =\left[ 
\begin{array}{cc}
a_{11} & 0 \\ 
0 & a_{22}%
\end{array}
\right] \left[ 
\begin{array}{c}
y_{1}\left( t\right) \\ 
y_{2}\left( t\right)%
\end{array}
\right]  \label{autonomuos}
\end{equation}
\begin{equation}
\rightarrow \left[ 
\begin{array}{c}
y_{1}\left( t\right) \\ 
y_{2}\left( t\right)%
\end{array}
\right] =\left[ 
\begin{array}{cc}
e^{a_{11}t} & 0 \\ 
0 & e^{a_{22}t}%
\end{array}
\right] \left[ 
\begin{array}{c}
c_{1} \\ 
c_{2}%
\end{array}
\right] .
\end{equation}

The system in (\ref{autonomuos}) has an important property, time has no
direct effect on the law-of-motion. Given knowledge of $y_{1}\left( t\right) 
$ and $y_{2}\left( t\right) $, $\dot{y}_{1}\left( t\right) $ and $\dot{y}%
_{2}\left( t\right) $ are fully determined, \emph{regardless of t. }A system
that has this property is called \emph{autonomous.} The system (\ref%
{diffsystem1}) is not an autonomous system, unless $\mathbf{p}\left(
t\right) $ is constant. Note that the homogeneous solution is always
autonomous, given, of course, that the parameters are constant.

The behavior of an autonomous system can be depicted in a graph, a phase
diagram. We see in the phase diagram that if the roots are stable, i.e.,
negative, the homogeneous part always goes to zero as $t$ goes to infinity.
With only one root stable, there is just one stable path.

\FRAME{itbpF}{6.1575in}{2.5002in}{0in}{}{}{slide6.emf}{\special{language
"Scientific Word";type "GRAPHIC";maintain-aspect-ratio TRUE;display
"USEDEF";valid_file "F";width 6.1575in;height 2.5002in;depth
0in;original-width 10.0016in;original-height 7.5135in;cropleft
"0.0377";croptop "1";cropright "0.9620";cropbottom "0.5037";filename
'Figures/Slide6.EMF';file-properties "XNPEU";}}

The fact that it is trivial to solve a diagonal system, suggests a way of
finding the solution to the homogeneous part of (\ref{diffsystem1}). Suppose
we can make a transformation of the variables so that the transformed system
is diagonal. Start by defining the new set of variables 
\begin{eqnarray}
\left[ 
\begin{array}{c}
x_{1}\left( t\right) \\ 
x_{2}\left( t\right)%
\end{array}
\right] &\equiv &\mathbf{B}\left[ 
\begin{array}{c}
y_{1}\left( t\right) \\ 
y_{2}\left( t\right)%
\end{array}
\right] \rightarrow \\
\left[ 
\begin{array}{c}
\dot{x}_{1}\left( t\right) \\ 
\dot{x}_{2}\left( t\right)%
\end{array}
\right] &\equiv &\mathbf{B}\left[ 
\begin{array}{c}
\dot{y}_{1}\left( t\right) \\ 
\dot{y}_{2}\left( t\right)%
\end{array}
\right] =\mathbf{BA}\left[ 
\begin{array}{c}
y_{1}\left( t\right) \\ 
y_{2}\left( t\right)%
\end{array}
\right] =\mathbf{BAB}^{-1}\left[ 
\begin{array}{c}
x_{1}\left( t\right) \\ 
x_{2}\left( t\right)%
\end{array}
\right]
\end{eqnarray}

If we can find a $\mathbf{B}$ such that $\mathbf{BAB}^{-1}$ is diagonal we
are half way. The solutions for $\mathbf{x}$ is then 
\begin{eqnarray}
&& \\
\left[ 
\begin{array}{c}
x_{1}\left( t\right) \\ 
x_{2}\left( t\right)%
\end{array}%
\right] &=&\left[ 
\begin{array}{cc}
e^{r_{1}t} & 0 \\ 
0 & e^{r_{2}t}%
\end{array}%
\right] \left[ 
\begin{array}{c}
c_{1} \\ 
c_{2}%
\end{array}%
\right]
\end{eqnarray}

where $r_{1,2}$ are the diagonal terms of the matrix $\mathbf{BAB}^{-1}$.The
solution for $\mathbf{y}$ then follows from the definition of $\mathbf{x}$%
\begin{equation}
\left[ 
\begin{array}{c}
y_{1}\left( t\right) \\ 
y_{2}\left( t\right)%
\end{array}
\right] =\mathbf{B}^{-1}\left[ 
\begin{array}{cc}
e^{r_{1}t} & 0 \\ 
0 & e^{r_{2}t}%
\end{array}
\right] \left[ 
\begin{array}{c}
c_{1} \\ 
c_{2}%
\end{array}
\right] .
\end{equation}

From linear algebra we know that $\mathbf{B}^{-1}$ is the matrix of
eigenvectors of $\mathbf{A}$ and that the diagonal terms of $\mathbf{BAB}%
^{-1}$ are the corresponding eigenvalues. The eigenvalues are given by the
characteristic equation of $\mathbf{A}$%
\begin{eqnarray}
\left\vert \left[ 
\begin{array}{cc}
a_{11}-r & a_{12} \\ 
a_{21} & a_{22}-r%
\end{array}%
\right] \right\vert &=&0 \\
&\rightarrow &\left( a_{11}-r\right) \left( a_{22}-r\right) -a_{12}a_{21}=0
\\
r^{2}-rTr\left( \mathbf{A}\right) +\left\vert \mathbf{A}\right\vert &=&0,
\end{eqnarray}%
where $Tr\left( \mathbf{A}\right) $ is the trace of $\mathbf{A.}$ The only
crux is that we need the roots to be distinct, otherwise $\mathbf{B}^{-1}$
is not always invertible. Distinct roots imply that $\mathbf{B}^{-1}$ is
invertible. (If $\mathbf{A}$ is symmetric $\mathbf{B}^{-1}$ is also
invertible.)

Let us draw a phase diagram with the eigenvectors of $\mathbf{A}$. The
dynamic system behaves as the diagonal one but the eigenvectors have
replaced the standard, orthogonal axes.

\FRAME{itbpF}{3.8043in}{2.5002in}{0in}{}{}{slide7.emf}{\special{language
"Scientific Word";type "GRAPHIC";maintain-aspect-ratio TRUE;display
"USEDEF";valid_file "F";width 3.8043in;height 2.5002in;depth
0in;original-width 10.2385in;original-height 7.5135in;cropleft "0";croptop
"1";cropright "0.5563";cropbottom "0.5037";filename
'Figures/Slide7.EMF';file-properties "XNPEU";}}

What is remaining is to find a particular solution of the complete system.
One way is here to use the method of undetermined coefficients. Sometimes,
we can find a steady state of the system, i.e., a point where the time
derivatives are all zero. This is easy if the exogenous part is constant. We
then set the differential equal to zero so 
\begin{eqnarray}
0 &=&\mathbf{Ay}\left( t\right) +\mathbf{p} \\
&\rightarrow &\mathbf{y}_{p}\left( t\right) =\mathbf{y}^{ss}=-\mathbf{A}^{-1}%
\mathbf{p.}
\end{eqnarray}

Given that we know the value of $\mathbf{y}\left( 0\right) $ we can now give
the general solution. The formula is given in matrix form and is valid for
any dimension of the system. First define

\begin{equation}
\mathbf{r}\left( t\right) \equiv \left[ 
\begin{array}{cccc}
e^{r_{1}t} & 0 & .. & 0 \\ 
0 & e^{r_{2}t} & .. & 0 \\ 
. & . & . & . \\ 
0 & 0 & . & e^{r_{n}t}%
\end{array}
\right] ,
\end{equation}
then we have 
\begin{eqnarray}
\mathbf{y}\left( t\right) &=&\mathbf{B}^{-1}\mathbf{x}\left( t\right) +%
\mathbf{y}^{ss}=\mathbf{B}^{-1}\mathbf{r}\left( t\right) \mathbf{c+y}^{ss} \\
\mathbf{y}\left( 0\right) &=&\mathbf{B}^{-1}\mathbf{x}\left( 0\right) +%
\mathbf{y}^{ss}=\mathbf{B}^{-1}\mathbf{c+y}^{ss} \\
\mathbf{c} &=&\mathbf{B}\left( \mathbf{\mathbf{y}\left( 0\right) -y}%
^{ss}\right) \\
&\rightarrow &\mathbf{y}\left( t\right) =\mathbf{B}^{-1}\mathbf{r}\left(
t\right) \mathbf{B}\left( \mathbf{\mathbf{y}\left( 0\right) -y}^{ss}\right) 
\mathbf{+y}^{ss}.
\end{eqnarray}

The method outlined above works also in the case of complex roots of the
characteristic equation. If the roots are $a\pm bi$ we have 
\begin{equation}
\mathbf{y}\left( t\right) =\mathbf{B}^{-1}\left[ 
\begin{array}{cc}
e^{\left( a+bi\right) t} & 0 \\ 
0 & e^{\left( a-bi\right) t}%
\end{array}
\right] \mathbf{c+y}^{ss}
\end{equation}
Example;

\begin{eqnarray}
\left[ 
\begin{array}{c}
\dot{y}_{1}\left( t\right) \\ 
\dot{y}_{2}\left( t\right)%
\end{array}
\right] &=&\left[ 
\begin{array}{cc}
-1 & -1 \\ 
1 & -1%
\end{array}
\right] \left[ 
\begin{array}{c}
y_{1}\left( t\right) \\ 
y_{2}\left( t\right)%
\end{array}
\right] +\left[ 
\begin{array}{c}
1 \\ 
1%
\end{array}
\right] \\
r_{1,2} &=&-1\pm i \\
\mathbf{B}^{-1} &=&\left[ 
\begin{array}{cc}
i & -i \\ 
1 & 1%
\end{array}
\right] .
\end{eqnarray}

So, 
\begin{gather}
\rightarrow \mathbf{y}\left( t\right) = \\
\left[ 
\begin{array}{cc}
i & -i \\ 
1 & 1%
\end{array}
\right] \left[ 
\begin{array}{cc}
e^{-t}\left( \cos t+i\sin t\right) & 0 \\ 
0 & e^{-t}\left( \cos t-i\sin t\right)%
\end{array}
\right] \left[ 
\begin{array}{c}
c_{1} \\ 
c_{2}%
\end{array}
\right] \mathbf{+y}^{ss} \\
=e^{-t}\left[ 
\begin{array}{c}
i\left( c_{1}-c_{2}\right) \cos t-\left( c_{1}+c_{2}\right) \sin t \\ 
\left( c_{1}+c_{2}\right) \cos t+i\left( c_{1}-c_{2}\right) \sin t%
\end{array}
\right] \mathbf{+y}^{ss} \\
=e^{-t}\left[ 
\begin{array}{c}
\tilde{c}_{1}\cos t+\tilde{c}_{2}\sin t \\ 
-\tilde{c}_{2}\cos t+\tilde{c}_{1}\sin t%
\end{array}
\right] \mathbf{+y}^{ss}
\end{gather}

The steady-state is found from, 
\begin{eqnarray}
\mathbf{0} &=&\left[ 
\begin{array}{cc}
-1 & -1 \\ 
1 & -1%
\end{array}
\right] \left[ 
\begin{array}{c}
y_{1}\left( t\right) \\ 
y_{2}\left( t\right)%
\end{array}
\right] +\left[ 
\begin{array}{c}
1 \\ 
1%
\end{array}
\right] \\
\left[ 
\begin{array}{c}
y_{1}^{ss} \\ 
y_{2}^{ss}%
\end{array}
\right] &=&-\left[ 
\begin{array}{cc}
-1 & -1 \\ 
1 & -1%
\end{array}
\right] ^{-1}\left[ 
\begin{array}{c}
1 \\ 
1%
\end{array}
\right] =\left[ 
\begin{array}{c}
0 \\ 
1%
\end{array}
\right] .
\end{eqnarray}

If we know $\mathbf{y}\left( 0\right) $, and using $\cos 0=1$ and $\sin 0=0,$%
we get 
\begin{gather}
\left[ 
\begin{array}{c}
y_{1}\left( 0\right) \\ 
y_{2}\left( 0\right)%
\end{array}
\right] =\left[ 
\begin{array}{c}
\tilde{c}_{1} \\ 
-\tilde{c}_{2}%
\end{array}
\right] \mathbf{+}\left[ 
\begin{array}{c}
0 \\ 
1%
\end{array}
\right] \\
\rightarrow \mathbf{y}\left( t\right) =e^{-t}\left[ 
\begin{array}{c}
y_{1}\left( 0\right) \cos t-\left( y_{2}\left( 0\right) -1\right) \sin t \\ 
\left( y_{2}\left( 0\right) -1\right) \cos t+y_{1}\left( 0\right) \sin t%
\end{array}
\right] \mathbf{+}\left[ 
\begin{array}{c}
0 \\ 
1%
\end{array}
\right]
\end{gather}

A phase-diagram of this system is an inward spiral.

\FRAME{itbpF}{3.2534in}{2.4984in}{0in}{}{}{slide8.emf}{\special{language
"Scientific Word";type "GRAPHIC";maintain-aspect-ratio TRUE;display
"USEDEF";valid_file "F";width 3.2534in;height 2.4984in;depth
0in;original-width 10.1866in;original-height 7.5299in;cropleft "0";croptop
"1";cropright "0.4786";cropbottom "0.5038";filename
'Figures/Slide8.EMF';file-properties "XNPEU";}}

How would it look like if the real part of the root was $-1$ or $+1$?

\subsubsection{Equivalent Systems}

In the repeated root case the matrix of the eigenvectors may be singular, so
that we cannot find $\mathbf{B}^{-1}$. Then we use the method of equivalent
systems, described in this section.

A linear \emph{nth} order differential equation is equivalent to a system of 
$n$ first order differential equations. Consider,\qquad \qquad 
\begin{equation}
\dddot{y}\left( t\right) +a_{1}\ddot{y}\left( t\right) +a_{2}\dot{y}\left(
t\right) +a_{3}y\left( t\right) =R\left( t\right)
\end{equation}

We can transform this into a system of first order differential equation by
defining 
\begin{eqnarray}
\dot{y}\left( t\right) &\equiv &x_{1}\left( t\right) , \\
\ddot{y}\left( t\right) &\equiv &x_{2}\left( t\right) =\dot{x}_{1}\left(
t\right) , \\
\dddot{y}\left( t\right) &=&\dot{x}_{2}\left( t\right) ,
\end{eqnarray}
giving 
\begin{equation}
\left[ 
\begin{array}{c}
\dot{y}\left( t\right) \\ 
\dot{x}_{1}\left( t\right) \\ 
\dot{x}_{2}\left( t\right)%
\end{array}
\right] =\left[ 
\begin{array}{ccc}
0 & 1 & 0 \\ 
0 & 0 & 1 \\ 
-a_{3} & -a_{2} & -a_{1}%
\end{array}
\right] \left[ 
\begin{array}{c}
y\left( t\right) \\ 
x_{1}\left( t\right) \\ 
x_{2}\left( t\right)%
\end{array}
\right] +\left[ 
\begin{array}{c}
0 \\ 
0 \\ 
R\left( t\right)%
\end{array}
\right]
\end{equation}

Since the equations are equivalent they consequently have the same
solutions. Sometimes one of the transformations is more convenient to solve.
Let us also transform a two dimensional system into a second order equation

\begin{equation}
\left[ 
\begin{array}{c}
\dot{x}_{1}\left( t\right) \\ 
\dot{x}_{2}\left( t\right)%
\end{array}
\right] =\left[ 
\begin{array}{cc}
a_{11} & a_{12} \\ 
a_{21} & a_{22}%
\end{array}
\right] \left[ 
\begin{array}{c}
x_{1}\left( t\right) \\ 
x_{2}\left( t\right)%
\end{array}
\right] +\left[ 
\begin{array}{c}
k_{1} \\ 
k_{2}%
\end{array}
\right] .  \label{simplesys}
\end{equation}

First use the first equation to express $x_{2}$ and then take the time
derivative of the first. Then we can eliminate $x_{2}$ and its time
derivative

\begin{gather}
\dot{x}_{1}\left( t\right) =a_{11}x_{1}\left( t\right) +a_{12}x_{2}\left(
t\right) +k_{1} \\
x_{2}\left( t\right) =\frac{\dot{x}_{1}\left( t\right) -a_{11}x_{1}\left(
t\right) -k_{1}}{a_{12}} \\
\dot{x}_{2}\left( t\right) =\frac{\ddot{x}_{1}\left( t\right) -a_{11}\dot{x}%
_{1}\left( t\right) }{a_{12}} \\
\frac{\ddot{x}_{1}\left( t\right) -a_{11}\dot{x}_{1}\left( t\right) }{a_{12}}%
=a_{21}x_{1}\left( t\right) +a_{22}\frac{\dot{x}_{1}\left( t\right)
-a_{11}x_{1}\left( t\right) -k_{1}}{a_{12}}+k_{2} \\
\ddot{x}_{1}\left( t\right) -\left( a_{11}+a_{22}\right) \dot{x}_{1}\left(
t\right) +\left( a_{22}a_{11}-a_{12}a_{21}\right) x_{1}\left( t\right)
=-a_{22}k_{1}+a_{12}k_{2}
\end{gather}

Note that the characteristic equation of this second order equation is the
same as the one for the system in (\ref{simplesys}). Consequently, the roots
and thus the dynamics are identical.

\subsection{Non-linear systems and Linearization\label{phasediagram}}

Phase diagrams are convenient to analyze the behavior of a 2 dimensional
system qualitatively that we cannot or prefer not to solve explicitly. E.g.,
if\qquad

\begin{eqnarray}
\dot{c}\left( t\right) &=&g_{1}\left( c\left( t\right) ,k\left( t\right)
\right) \\
\dot{k}\left( t\right) &=&g_{2}\left( c\left( t\right) ,k\left( t\right)
\right)
\end{eqnarray}

The first step here is to find the two curves in the $c,k$-space where $c$
and $k$, respectively are constant. Setting the time derivatives equal to
zero defines two relations between $c$ and $k$, which we denote by $G_{1}$
and $G_{2}.$%
\begin{eqnarray}
g_{1}\left( c\left( t\right) ,k\left( t\right) \right) &=&0\rightarrow
c=G_{1}\left( k\right) \\
g_{2}\left( c\left( t\right) ,k\left( t\right) \right) &=&0\rightarrow
c=G_{2}\left( k\right)
\end{eqnarray}%
We then draw these curves in the $c,k$-space. For example, you will in the
macro course analyze the Ramsey optimal consumption problem, where output is 
$f\left( k\right) ,$ interest rate is $f^{\prime }\left( k\right) ,$ the
subjective discount rate is $\theta $ and utility is $u\left( c\right) $
where $u$ and $f$ are assumed to be concave functions. The model will
deliver the following system of differential equations 
\begin{eqnarray}
\dot{c}\left( t\right) &=&-\frac{u^{\prime }\left( c\left( t\right) \right) 
}{u^{\prime \prime }\left( c\left( t\right) \right) }\left( f^{\prime
}\left( k\left( t\right) \right) -\theta \right)  \label{Ramsey} \\
\dot{k}\left( t\right) &=&f\left( k\left( t\right) \right) -c\left( t\right)
.
\end{eqnarray}%
Setting the time derivatives to zero we get \ \qquad 
\begin{eqnarray}
f^{\prime }\left( k\left( t\right) \right) &=&\theta , \\
c &=&f\left( k\left( t\right) \right) .
\end{eqnarray}

Draw these curves in the $c,k$ space

\FRAME{ftbpF}{4.651in}{3.2041in}{0pt}{}{}{slide10.emf}{\special{language
"Scientific Word";type "GRAPHIC";maintain-aspect-ratio TRUE;display
"USEDEF";valid_file "F";width 4.651in;height 3.2041in;depth
0pt;original-width 10.0016in;original-height 7.5135in;cropleft "0";croptop
"1";cropright "0.6971";cropbottom "0.3626";filename
'Figures/Slide10.EMF';file-properties "XNPEU";}}

We then have to find the signs of time derivatives, and above and below
their respective zero motion curves. From (\ref{Ramsey}), we see that 
\begin{eqnarray}
\frac{\partial \dot{c}\left( t\right) }{\partial k} &=&-\frac{u^{\prime
}\left( c\right) }{u^{\prime \prime }\left( c\right) }f^{\prime \prime
}\left( k\right) <0 \\
\frac{\partial \dot{k}\left( t\right) }{\partial c} &=&-1.
\end{eqnarray}

This means that $\dot{c}$ is positive to the left of and negative to the
right of the curve $\dot{c}=0$. For $\dot{k}$, we find that it is positive
below and negative above $\dot{k}=0$. Then draw these motions as arrows in
the phase diagram. Note that no paths ever can cross.

\FRAME{itbpF}{4.651in}{3.0026in}{0in}{}{}{slide11.emf}{\special{language
"Scientific Word";type "GRAPHIC";maintain-aspect-ratio TRUE;display
"USEDEF";valid_file "F";width 4.651in;height 3.0026in;depth
0in;original-width 10.0016in;original-height 7.5135in;cropleft "0";croptop
"1";cropright "0.6971";cropbottom "0.4029";filename
'Figures/Slide11.EMF';file-properties "XNPEU";}}

We conclude that this system has saddle point characteristics and thus has
only one stable trajectory towards the steady state.

The behavior close to the steady state should also be evaluated by means of
linearization around the steady state. We do that by approximating in the
following way 
\begin{equation}
\left[ 
\begin{array}{c}
\dot{c}\left( t\right) \\ 
\dot{k}\left( t\right)%
\end{array}%
\right] \approx \left[ 
\begin{array}{cc}
\frac{\partial g_{1}\left( c^{ss},k^{ss}\right) }{\partial c} & \frac{%
\partial g_{1}\left( c^{ss},k^{ss}\right) }{\partial k} \\ 
\frac{\partial g_{2}\left( c^{ss},k^{ss}\right) }{\partial c} & \frac{%
\partial g_{2}\left( c^{ss},k^{ss}\right) }{\partial c}%
\end{array}%
\right] \left[ 
\begin{array}{c}
c\left( t\right) -c^{ss} \\ 
k\left( t\right) -k^{ss}%
\end{array}%
\right] .
\end{equation}%
with an obvious generalization to higher dimensions.

We now evaluate the roots of the matrix of derivatives. In the example we
find that the coefficient matrix is

\begin{eqnarray}
&&\left[ 
\begin{array}{cc}
\frac{\partial g_{1}\left( c^{ss},k^{ss}\right) }{\partial c} & \frac{%
\partial g_{1}\left( c^{ss},k^{ss}\right) }{\partial k} \\ 
\frac{\partial g_{2}\left( c^{ss},k^{ss}\right) }{\partial c} & \frac{%
\partial g_{2}\left( c^{ss},k^{ss}\right) }{\partial c}%
\end{array}
\right] \\
&=&\left[ 
\begin{array}{cc}
-\left( f^{\prime }-\theta \right) \frac{\partial }{\partial c}\left( \frac{%
u^{\prime }}{u^{\prime \prime }}\right) & -\frac{u^{\prime }}{u^{\prime
\prime }}f^{\prime \prime } \\ 
-1 & f^{\prime }%
\end{array}
\right] ,
\end{eqnarray}
where all functions are evaluated at the steady state. There, $f^{\prime
}=\theta $, implying that the matrix simplifies to 
\begin{equation}
\left[ 
\begin{array}{cc}
0 & -\frac{u^{\prime }}{u^{\prime \prime }}f^{\prime \prime } \\ 
-1 & f^{\prime }%
\end{array}
\right]
\end{equation}
with eigenvalues 
\begin{equation}
\frac{1}{2}\left( f^{\prime }\pm \sqrt{(f^{\prime })^{2}+4\frac{u^{\prime }}{%
u^{\prime \prime }}f^{\prime \prime }}\right)
\end{equation}
which clearly are of opposite signs.

\subsection{Example: Steady-state asset distributions}

In this example, we use our derived skills to find the steady-state wealth
distribution in a simple model. As we will see, here the model gives us a
system of differential equations for the wealth distribution that we can
solve easily. We will see that the methods we have learned work as well when
the differential equations are in wealth, rather than time.

In (Hassler and Mora, JPub, 1999), we analyze preferences over unemployment
insurance in a very simple continuous time economy where agents can borrow
and save at rate $r$. In the model, employed individuals earn a wage $w$ per
unit of time and become unemployed with an instantaneous probability $q.$
This means that over a small (infinitesimal) interval of time $dt,$ the
probability of becoming unemployed is $qdt.$ Unemployed individuals receive
a flow of unemployment benefits $b$ and become rehired with instantaneous
probability $hdt.$ In addition, we assume that there is a instantaneous
death probability of $\delta $ and that there is an inflow of newborn
unemployed with zero assets of $\delta $ so that the total population is
constant.

In the paper, we show that if individuals have CARA utility ($U=-e^{\gamma
c} $) and wages and benefits are constant, individuals choose constant
savings amounts for each of the two employment states. Denoting these $s_{e}$
and $s_{u}$, where $s_{e}>0$ and $s_{u}<0$, we have that individual asset
accumulation for the two types, conditional on surviving is 
\begin{equation}
\dot{A}_{t}=s_{e},
\end{equation}%
for employed and 
\begin{equation}
\dot{A}_{t}=s_{u}
\end{equation}%
for unemployed.

In the paper, we don't calculate the steady state wealth distribution of
assets. That's the purpose of this exercise. First, we calculate the steady
state share of unemployed, $\mu _{u}$. For this purpose, we note that in
steady state, the inflow and the outflow to the stock of unemployed must be
constant. Thus,

\begin{eqnarray}
\left( 1-\mu _{u}\right) q+\delta &=&\mu _{u}\left( h+\delta \right) \\
&\rightarrow &\mu _{u}=\frac{q+\delta }{q+h+\delta }
\end{eqnarray}

Now, we want to calculate the steady state distribution of assets among
employed and unemployed in this economy. Let us denote these densities, by $%
f_{e}(A)$ and $f_{u}(A).$ We will derive these by solving a system of
differential equations. Consider first the number (density) of employed
individuals with assets $A\neq 0$, given by $f_{e}\left( A\right) $. Over a
small period $dt$, a number $f_{e}(A)\left( 1-\left( \delta +q\right)
dt\right) $ of them remain employed and alive. Over the same time, a number $%
f_{u}(A)hdt$ of unemployed with assets $A$ become hired. Finally, over the
period $dt$ these individuals add to their assets an amount $s_{e}dt.$
Writing this down yields, 
\begin{eqnarray}
f_{e}(A+s_{e}dt) &=&f_{e}(A)\left( 1-\left( \delta +q\right) dt\right)
+f_{u}(A)hdt, \\
f_{e}(A)+f_{e}^{\prime }(A)s_{e}dt &=&f_{e}(A)\left( 1-\left( \delta
+q\right) dt\right) +f_{u}(A)hdt \\
f_{e}^{\prime }(A) &=&-f_{e}(A)\frac{\delta +q}{s_{e}}+f_{u}(A)\frac{h}{s_{e}%
}.
\end{eqnarray}%
By the same reasoning,

\begin{eqnarray}
f_{u}(A+s_{e}dt) &=&f_{u}(A)\left( 1-\left( \delta +h\right) dt\right)
+f_{e}(A)qdt, \\
f_{u}(A)+f_{u}^{\prime }(A)s_{e}dt &=&f_{u}(A)\left( 1-\left( \delta
+h\right) dt\right) +f_{e}(A)qdt \\
f_{u}^{\prime }(A) &=&-f_{u}(A)\frac{\delta +h}{s_{u}}+f_{e}(A)\frac{q}{s_{u}%
},
\end{eqnarray}%
yielding the system 
\begin{eqnarray}
\left[ 
\begin{array}{c}
f_{e}^{\prime }(A) \\ 
f_{u}^{\prime }(A)%
\end{array}%
\right] &=&\left[ 
\begin{array}{cc}
-\frac{\delta +q}{s_{e}} & \frac{h}{s_{e}} \\ 
\frac{q}{s_{u}} & -\frac{\delta +h}{s_{u}}%
\end{array}%
\right] \left[ 
\begin{array}{c}
f_{e}(A) \\ 
f_{u}(A)%
\end{array}%
\right] . \\
&\equiv &\mathbf{A}\left[ 
\begin{array}{c}
f_{e}(A) \\ 
f_{u}(A)%
\end{array}%
\right]
\end{eqnarray}

Let us now, assign some numbers to the parameters. Say $\delta
=1/40,h=2,q=1/5,s_{e}=1$ and $s_{u}=-2$, all measured as probabilities per
year. Then, we can calculate the roots and the eigenvalues 
\begin{eqnarray}
r_{1} &=&\frac{63}{160}+\frac{1}{160}\sqrt{4681}\approx .821 \\
r_{2} &=&\frac{63}{160}-\frac{1}{160}\sqrt{4681}\approx -0.0339 \\
\mathbf{B}^{-1} &=&\left[ 
\begin{array}{cc}
1 & 1 \\ 
\frac{99}{320}+\frac{1}{320}\sqrt{4681} & \frac{99}{320}-\frac{1}{320}\sqrt{%
4681}%
\end{array}%
\right] \approx \left[ 
\begin{array}{cc}
1.0 & 1.0 \\ 
.523 & 0.0956%
\end{array}%
\right] .
\end{eqnarray}

Thus, except at $A=0$, where there is an inflow of newborn, that we have not
considered, we have

\begin{gather}
\left[ 
\begin{array}{c}
f_{e}(A) \\ 
f_{u}(A)%
\end{array}
\right] =\mathbf{B}^{-1}\mathbf{r}\left( t\right) \mathbf{c} \\
=\left[ 
\begin{array}{cc}
1.0 & 1.0 \\ 
.523 & 0.0956%
\end{array}
\right] \left[ 
\begin{array}{cc}
e^{.821A} & 0 \\ 
0 & e^{-0.0339A}%
\end{array}
\right] \left[ 
\begin{array}{c}
c_{1} \\ 
c_{2}%
\end{array}
\right] \\
=\left[ 
\begin{array}{c}
c_{1}e^{.821A}+c_{2}e^{-.0339A} \\ 
.523c_{1}e^{.821A}+.0956c_{2}e^{-.0339A}%
\end{array}
\right]
\end{gather}

Now, we first note that $f_{e}(A)$ and $f_{u}(A)$ cannot be explosive in any
of the directions. That would violate that these functions are densities,
i.e., the sum of their respective integrals over the real line must be
unity. Thus, for $A<0$, $c_{2}$ must be zero and for $A>0,$ $c_{1}=0.$
Furthermore, we know 
\begin{eqnarray}
\int_{-\infty }^{\infty }f_{e}(A)dA &=&1-\mu _{u}=\frac{h}{q+h+\delta }%
\approx 0.899 \\
\int_{-\infty }^{\infty }f_{u}(A)dA &=&\mu _{u}=\frac{q+\delta }{q+h+\delta }%
\approx 0.101.
\end{eqnarray}

Using this, we can calculate the integrations constants. 
\begin{eqnarray}
\int_{-\infty }^{\infty }f_{e}(A)dA &=&\int_{-\infty
}^{0}c_{1}e^{.821A}dA+\int_{0}^{\infty }c_{2}e^{-.0339A}dA \\
&=&\frac{c_{1}}{0.821}+\frac{c_{2}}{.0339}=0.899 \\
\int_{-\infty }^{\infty }f_{u}(A)dA &=&\int_{-\infty
}^{0}.523c_{1}e^{.821A}dA+\int_{0}^{\infty }.0956e^{-.0\,\allowbreak 339A}dA
\\
&=&\frac{.523c_{1}}{0.821}+\frac{0.0956c_{2}}{.0339}=0.101
\end{eqnarray}

Yielding,%
\begin{eqnarray}
c_{1} &=&0.0289, \\
c_{2} &=&0.0293.
\end{eqnarray}

This concludes our calculations,

\begin{eqnarray}
f_{e}(A) &=&\left\{ 
\begin{array}{c}
0.0289e^{.\,821A}\text{ for }A<0 \\ 
0.0293e^{-.0339A}\text{ for }A>0%
\end{array}%
\right. , \\
f_{u}(A) &=&\left\{ 
\begin{array}{c}
0.0151e^{.821A}\text{ for }A<0 \\ 
0.00280e^{-.0\,339A}\text{ for }A>0%
\end{array}%
\right. .
\end{eqnarray}%
\FRAME{itbpFX}{2.9996in}{2.3896in}{0in}{}{}{Plot}{\special{language
"Scientific Word";type "MAPLEPLOT";width 2.9996in;height 2.3896in;depth
0in;display "USEDEF";plot_snapshots TRUE;mustRecompute FALSE;lastEngine
"Maple";xmin "-5.000000";xmax "0";xviewmin "-5.000000";xviewmax "0";yviewmin
"0";yviewmax "0.04";viewset"XY";plottype 4;numpoints 49;plotstyle
"patch";axesstyle "normal";xis \TEXUX{A};var1name \TEXUX{$A$};function
\TEXUX{$0.0289e^{.\,\allowbreak 821A}$};linecolor "black";linestyle
1;pointstyle "point";linethickness 1;lineAttributes "Solid";var1range
"-5.000000,0";num-x-gridlines 49;curveColor
"[flat::RGB:0000000000]";curveStyle "Line";function
\TEXUX{$0.0151e^{.\,\allowbreak 821A}$};linecolor "black";linestyle
3;pointstyle "point";linethickness 1;lineAttributes "Dots";var1range
"-5,0";num-x-gridlines 49;curveColor "[flat::RGB:0000000000]";curveStyle
"Line";rangeset"X";valid_file "T";tempfilename
'KBKCSE02.wmf';tempfile-properties "XPR";}}\FRAME{itbpFX}{2.9996in}{2.3896in%
}{0in}{}{}{Plot}{\special{language "Scientific Word";type "MAPLEPLOT";width
2.9996in;height 2.3896in;depth 0in;display "USEDEF";plot_snapshots
TRUE;mustRecompute FALSE;lastEngine "Maple";xmin "0.000000";xmax
"100.000000";xviewmin "0.000000";xviewmax "100.000000";yviewmin "0";yviewmax
"0.04";viewset"XY";plottype 4;numpoints 49;plotstyle "patch";axesstyle
"normal";xis \TEXUX{A};var1name \TEXUX{$A$};function \TEXUX{$0.02\allowbreak
93e^{-.0\,\allowbreak 339A}$};linecolor "black";linestyle 1;pointstyle
"point";linethickness 1;lineAttributes "Solid";var1range
"0.000000,100.000000";num-x-gridlines 49;curveColor
"[flat::RGB:0000000000]";curveStyle "Line";function \TEXUX{$0.002\allowbreak
80e^{-.0\,\allowbreak 339A}$};linecolor "cyan";linestyle 3;pointstyle
"point";linethickness 1;lineAttributes "Dots";var1range
"0,100";num-x-gridlines 49;curveColor "[flat::RGB:0x00008080]";curveStyle
"Line";rangeset"X";valid_file "T";tempfilename
'KBKCSF03.wmf';tempfile-properties "XPR";}}

Plotting the densities, using a solid (dotted) line to denote wealth of
employed (unemployed), we see that unemployed are more concentrated to the
left. We can also calculate average assets among the two types from 
\begin{gather}
A_{e}=\int_{-\infty }^{\infty }Af_{e}(A)dA= \\
\int_{-\infty }^{0}A\cdot 0.0289e^{.821A}dA+\int_{0}^{\infty }A\cdot
0.0293e^{-.0\,\allowbreak 339A}dA\approx 25.45 \\
A_{u}=\int_{-\infty }^{\infty }Af_{u}(A)dA= \\
\int_{-\infty }^{0}A\cdot 0.0151e^{.821A}dA+\int_{0}^{\infty }A\cdot
0.00280e^{-.0\,339A}dA\approx 2.41
\end{gather}
\newpage

\section{Difference equations}

\subsection{Sums, forward and backward solutions}

\subsubsection{Sums vs. Integrals}

%TCIMACRO{\TeXButton{TeX field}{\setcounter{equation}{0}}}%
%BeginExpansion
\setcounter{equation}{0}%
%EndExpansion
Difference equations can be solved in ways very similar to how we solve
differential equations. First, we will look at the analogy to integrals.

Consider the difference equation 
\begin{equation}
A_{t}-A_{t-1}\equiv \Delta A_{t}=q_{t}  \label{differenceeq1}
\end{equation}%
where $\Delta A_{t}$ is the change in $A$ per unit (interval) of time. We
sum both sides from some date $t_{0}$ until $t$ to get

\begin{equation}
\sum_{s=t_{0}}^{t}\Delta A_{s}=A_{t}-A_{t_{0}-1}=\sum_{s=t_{0}}^{t}q_{s}.
\end{equation}

We can then write the solution as 
\begin{eqnarray}
\Delta A_{t} &=&q_{t} \\
&\rightarrow &A_{t}=\sum_{s=t_{0}}^{t}q_{s}+A_{t_{0}-1}.  \label{difference1}
\end{eqnarray}

As we see, this is very much like integrals 
\begin{eqnarray}
\frac{dA\left( t\right) }{dt} &=&q\left( t\right) \\
&\rightarrow &A\left( t\right) =\int_{t_{0}}^{t}q\left( s\right)
ds+A_{t_{0}},
\end{eqnarray}
and the relation between $q_{t}$ and $\sum_{s=t_{0}}^{t}q_{s}$ is the same
as between $q\left( t\right) $ and its primitive, since 
\begin{equation}
\Delta \sum_{s=t_{0}}^{t}q_{s}=q_{t}.
\end{equation}

\subsubsection{Forward and backward solutions}

The part $\sum_{s=t_{0}}^{t}q_{t}$ of the RHS of (\ref{difference1}) is
exogenous, i.e., independent of $A_{t}$. Sometimes, it converges in one or
both directions. Then we can write the solutions in another way. Suppose the
following limit exists

\begin{equation}
\lim_{T\rightarrow \infty }\sum_{s=-T}^{t}q_{s}\equiv \sum_{s=-\infty
}^{t}q_{s}.
\end{equation}
Then, the other part of RHS of (\ref{difference1}) should also have a
well-defined limit 
\begin{equation}
\lim_{T\rightarrow \infty }A_{-T}\equiv \underline{A},
\end{equation}
so that the solution is 
\begin{equation}
A_{t}=\sum_{s=-\infty }^{t}q_{s}+\underline{A}.  \label{backsol}
\end{equation}

Clearly, (\ref{backsol}) solves (\ref{differenceeq1}),

\begin{equation}
A_{t}-A_{t-1}=\sum_{s=-\infty }^{t}q_{s}+\underline{A}-\sum_{s=-\infty
}^{t-1}q_{s}-\underline{A}=q_{t}.
\end{equation}%
If the solution in (\ref{backsol}) exists, it is called the \emph{backward
solution. }

Analogously, the limit 
\begin{equation}
\lim_{T\rightarrow \infty }\sum_{s=t}^{T}q_{s}\equiv \sum_{s=t}^{\infty
}q_{s},
\end{equation}%
might exist, in which case

\begin{equation}
\lim_{T\rightarrow \infty }A_{T}\equiv \bar{A}
\end{equation}
also exists. Then, we have 
\begin{eqnarray}
\lim_{T\rightarrow \infty }\left( A_{T}-A_{t}\right) &=&\sum_{s=t+1}^{\infty
}q_{s}, \\
&\rightarrow &A_{t}=\bar{A}-\sum_{s=t+1}^{\infty }q_{s},
\end{eqnarray}
which is called the \emph{forward solution}.

Example. Suppose $q_{t}$ follows a simple $AR(1)$ process

\begin{eqnarray}
\Delta A_{t} &=&q_{t}, \\
q_{t} &=&rq_{t-1.}
\end{eqnarray}

If $r<1,$ we can use the forward solution and if $r>1$, the backward
solution works. In the former case, 
\begin{eqnarray}
\sum_{s=t+1}^{\infty }q_{s} &=&\sum_{s=0}^{\infty }q_{t+1}r^{s}=\frac{q_{t+1}%
}{1-r}, \\
&\rightarrow &A_{t}=\bar{A}-\frac{q_{t+1}}{1-r},
\end{eqnarray}%
where $\bar{A}$ is determined from, for example, an initial condition. In
the latter case, 
\begin{eqnarray}
\sum_{s=-\infty }^{t}q_{s} &=&\sum_{s=0}^{\infty }q_{t}r^{-s}=q_{t}\frac{r}{%
r-1} \\
A_{t} &=&q_{t}\frac{r}{r-1}+\underline{A}.
\end{eqnarray}

\subsubsection{First order difference equations with constant coefficients}

A first order difference equation with constant coefficients has the
following form.

\begin{equation}
x_{t}-ax_{t-1}=c  \label{differenceeq2}
\end{equation}

As we see, the LHS is not a pure difference, as in (\ref{differenceeq1}), so
we cannot simply sum over $t$. Instead we rely on the following result.

\begin{result}
The general solution (the complete set of solutions) to a difference
equation is the general solution to the homogeneous part plus any particular
solution to the complete equation.

The general solution to the homogeneous first order difference equation with
coefficient $a$ can be written

\begin{equation}
x_{t}^{h}=Aa^{t}
\end{equation}

where $A$ is an arbitrary constant.
\end{result}

A particular solution is sometimes the steady state which exists for (\ref%
{differenceeq2}) if $a\neq 1.$%
\begin{eqnarray}
x^{ss}-ax^{ss} &=&c \\
&\rightarrow &x^{ss}=\frac{c}{1-a}.
\end{eqnarray}

\subsubsection{Non-constant RHS}

Now look at 
\begin{equation}
x_{t}-ax_{t-1}=q_{t}.  \label{differenceeq3}
\end{equation}

A way to solve (\ref{differenceeq3}) is to use $a^{-t}$ as the analogue to
the integrating factor. Multiplying both sides by $a^{-t}$ yields 
\begin{equation}
a^{-t}\left( x_{t}-ax_{t-1}\right) =a^{-t}q_{t}.
\end{equation}%
Now, we see that the LHS\ can be written $\Delta \left( a^{-t}x_{t}\right) $%
, implying 
\begin{eqnarray}
\Delta \left( a^{-t}x_{t}\right) &=&a^{-t}q_{t}, \\
&\Rightarrow &a^{-t}x_{t}=\sum_{s=t_{0}}^{t}a^{-s}q_{s}+A, \\
a^{-t}x_{t} &=&\sum_{s=t_{0}}^{t}a^{-s}q_{s}+A, \\
x_{t} &=&Aa^{t}+\sum_{s=t_{0}}^{t}a^{t-s}q_{s}.
\end{eqnarray}%
Again, we should verify that this satisfies our original difference equation.

\subsubsection{Stable growth -- the Solow growth model}

Often in macro, the variables in the model grow in a way that precludes the
existence of a steady state. However, some transformation of the variables
might possess a steady state. The simplest example of this is the Solow
growth model \ We will first solve the model in trending variables and then
do it in transformed (detrended) variables. As we know, the savings rate is
exogenous and denoted $S$ and the labor supply $N_{t}$ follows 
\begin{equation}
N_{t}=e^{g}N_{t-1}\approx \left( 1+g\right) N_{t-1}.
\end{equation}

There is one good, used for consumption and as capital, which follows the
law-of-motion 
\begin{equation}
K_{t+1}=Sf\left( K_{t},N_{t}\right) ,
\end{equation}
where $f\left( K_{t},N_{t}\right) $ is a concave production function. Let us
specify production as the CRS Cobb-Douglas function, 
\begin{eqnarray}
f\left( K_{t},N_{t}\right) &=&N_{t}^{1-\alpha }K_{t}^{\alpha },\text{ }
\label{Solow} \\
1 &>&\alpha >0.
\end{eqnarray}

Letting lower case letters denote natural logarithms, we have 
\begin{eqnarray}
k_{t+1} &=&s+\left( 1-\alpha \right) n_{t}+\alpha k_{t}, \\
k_{t+1}-\alpha k_{t} &=&s+\left( 1-\alpha \right) n_{t},
\end{eqnarray}%
and 
\begin{eqnarray}
\Delta n_{t} &=&g \\
&\rightarrow &n_{t}=\sum_{0}^{t}g+n=tg+n,
\end{eqnarray}%
for some constant $n.$

Here, we can guess on a particular solution of the same form as the RHS,
i.e., 
\begin{eqnarray}
k_{t} &=&k+tg_{k} \\
k+\left( t+1\right) g_{k}-\alpha \left( k+tg_{k}\right) &=&s+\left( 1-\alpha
\right) n_{t} \\
k\left( 1-\alpha \right) +t\left( 1-\alpha \right) g_{k}+g_{k} &=&s+\left(
1-\alpha \right) n+t\left( 1-\alpha \right) g
\end{eqnarray}

For this to hold for all $t$, we need 
\begin{eqnarray}
g_{k} &=&g, \\
k &=&\frac{s-g}{1-\alpha }+n.
\end{eqnarray}
giving a particular solution 
\begin{equation}
k_{p,t}=\frac{s-g}{1-\alpha }+n+tg.
\end{equation}
This solution is in economics often called a \emph{balanced growth path}, in
which all variables grow at the same rate. The complete solution is 
\begin{eqnarray}
k_{t} &=&Aa^{t}+\frac{s-g}{1-\alpha }+n+tg,  \label{balancedgrwth} \\
k_{t} &=&\left( k_{0}-\left( \frac{s-g}{1-\alpha }+n\right) \right) a^{t}+%
\frac{s-g}{1-\alpha }+n+tg.
\end{eqnarray}

In this case, a convenient alternative is to define a new variable, capital
per capita, $C_{t}\equiv K_{t}/N_{t}$ $\rightarrow $ $c_{t}=k_{t}-n_{t}$.
Using this, and dividing (\ref{Solow}) by $N_{t}$, we get 
\begin{eqnarray}
\frac{K_{t+1}}{N_{t}} &=&\frac{Sf\left( K_{t},N_{t}\right) }{N_{t}} \\
\frac{K_{t+1}N_{t+1}}{N_{t+1}N_{t}} &=&\frac{SN_{t}^{1-\alpha }K_{t}^{\alpha
}}{N_{t}}= \\
C_{t+1}e^{g} &=&S\left( \frac{K_{t}}{N_{t}}\right) ^{\alpha } \\
c_{t+1}+g &=&s+\alpha c_{t} \\
c_{t+1}-\alpha c_{t} &=&s-g \\
&\rightarrow &c^{ss}=\frac{s-g}{1-\alpha } \\
c_{t} &=&\left( c_{0}-\frac{s}{1-\alpha }\right) \alpha ^{t}+\frac{s}{%
1-\alpha },
\end{eqnarray}%
coinciding with the solution in (\ref{balancedgrwth}) but now expressed as
steady state in capital per capita, rather than a balanced growth path for
capital.

We can also look at the log of output per capita, 
\begin{eqnarray}
y_{t} &=&\alpha c_{t}, \\
y_{t} &=&\left( y_{0}-\alpha \frac{s}{1-\alpha }\right) \alpha ^{t}+\alpha 
\frac{s}{1-\alpha }.
\end{eqnarray}

This is the basis for the so-called growth regressions, pioneered by Barro, 
\begin{equation}
y_{t}-y_{0}=-y_{0}\left( 1-\alpha ^{t}\right) +\frac{s}{1-\alpha }\alpha
\left( 1-\alpha ^{t}\right) ,
\end{equation}
where growth over a sample period $0$ through $t$ is seen to depend
negatively on initial output and positively on savings. Also, 
\begin{eqnarray}
y_{t}-y_{t-1} &=&\left( y_{0}-\alpha \frac{s}{1-\alpha }\right) \alpha
^{t}+\alpha \frac{s}{1-\alpha } \\
&&-\left( y_{0}-\alpha \frac{s}{1-\alpha }\right) \alpha ^{t-1}-\alpha \frac{%
s}{1-\alpha } \\
&=&-\left( y_{0}-\alpha \frac{s}{1-\alpha }\right) \alpha ^{t-1}\left(
1-\alpha \right) \\
&=&\left( \alpha \frac{s}{1-\alpha }-y_{t-1}\right) \left( 1-\alpha \right) .
\end{eqnarray}

As we see, the growth rate (the log difference) of output per capita is a
fraction $\left( 1-\alpha \right) $ of the difference between the steady
state and current output per capita. Note that convergence is \emph{slower }%
the larger is the capital's share of output. In the limit when $\alpha
\rightarrow 1$, the model shows no convergence and has become an \emph{%
endogenous growth} model, with 
\begin{eqnarray}
\Delta k_{t} &=&s, \\
\Delta c_{t} &=&\Delta y_{t}=s-g.
\end{eqnarray}

\subsection{Linear difference equations of higher order}

\subsubsection{Higher order homogeneous difference equations with constant
coefficients}

Consider the homogeneous difference equation 
\begin{equation}
x_{t+n}+a_{1}x_{t+n-1}+...a_{n}x_{t}=0,  \label{diffpoly}
\end{equation}

\begin{definition}
The forward operator E is defined by 
\begin{equation}
E^{s}x_{t}\equiv Ex_{t+s}
\end{equation}
where $s$ is any integer, positive or negative.
\end{definition}

We can then write (\ref{diffpoly}) in a condensed polynomial form

\begin{equation}
P\left( E\right) x_{t}=0.
\end{equation}

We then have to find the roots of the equation

\begin{equation}
P\left( r\right) =r^{n}+a_{1}r^{n-1}+...a_{n}=0
\end{equation}

Each root contributes to the general solution with one term that is
independent of the others, exactly as with differential equations.

\begin{result}
Let $r_{s}$ denote the roots to the polynomial $P\left( r\right) $, i.e.,
all solutions to $P(r)=0$. Let the first $k\geq 0$ roots be distinct and the
remaining $l=n-k$ roots repeated. Then, the general solution to (\ref%
{diffpoly}) is 
\begin{equation}
x_{t}=c_{1}r_{1}^{t}+...+c_{k}r_{k}^{t}+r_{k+1}^{t}\left(
c_{k+1}+tc_{k+2}+...+t^{l-1}c_{k+l}\right)  \label{solutiondiffpoly}
\end{equation}
If there is more than one set of repeated roots, each set of size $m\geq 2$
contributes with $m$ linearly independent terms 
\begin{equation}
r_{k+1}^{t}\left( c_{k+1}+tc_{k+2}+...+t^{m-1}c_{k+m}\right) .
\end{equation}
\end{result}

In the case of complex roots, express the complex number in polar form%
\footnote{%
Note that the formula $\theta =\tan ^{-1}\frac{b}{a}$ is valid only for $%
\theta \in \left( -\pi /2,\pi /2\right] .$ For example, if $r_{1,2}=-\frac{1%
}{2}\pm \frac{1}{2}\sqrt{3}i$ we note that 
\begin{eqnarray*}
-\frac{1}{2}+\frac{1}{2}\sqrt{3}i &=&e^{-i\frac{4}{3}\pi }, \\
-\frac{1}{2}-\frac{1}{2}\sqrt{3}i &=&e^{i\frac{4}{3}\pi }.
\end{eqnarray*}%
} 
\begin{eqnarray}
r &=&a+bi=\left\vert r\right\vert \left( \cos \theta +i\sin \theta \right) ,
\\
\left\vert r\right\vert &=&\sqrt{a^{2}+b^{2},} \\
\theta &=&\tan ^{-1}\frac{b}{a}.
\end{eqnarray}

We then use the fact that 
\begin{eqnarray}
e^{i\theta } &=&\left( \cos \theta +i\sin \theta \right) \\
r &=&\left| r\right| e^{i\theta }\rightarrow r^{t}=\left| r\right|
^{t}e^{i\theta t}=\left| r\right| ^{t}\left( \cos t\theta +i\sin t\theta
\right)
\end{eqnarray}
to get for the complex conjugates $r_{1,2}=a\pm bi,$ 
\begin{equation}
c_{1}r_{1}^{t}+c_{2}r_{2}^{t}=\left| r\right| ^{t}\left( \tilde{c}_{1}\cos
t\theta +\tilde{c}_{2}\sin t\theta \right)  \label{solutioncomplex}
\end{equation}
which we can see is a generalization of (\ref{solutiondiffpoly}). Complex
roots thus give us oscillating solutions.

\subsubsection{Stability}

From the solutions (\ref{solutiondiffpoly}) and (\ref{solutioncomplex}), it
is clear that roots such that $\left| r\right| <1,$give converging terms.

\subsubsection{Higher order non-homogeneous difference equations with
constant coefficients}

\begin{result}
The general solution to a non-homogeneous difference equation with constant
coefficients is given by the general solution to the homogeneous part plus
any solution to the full equation.
\end{result}

To solve non-homogeneous equations we thus have to find particular solution
to the complete equation to add to the general solution of the homogeneous
part. The simplest non-homogeneous difference equation is 
\begin{equation}
P\left( x_{t}\right) =c.
\end{equation}

Here, we try a steady state 
\begin{eqnarray}
P\left( x^{ss}\right) &=&P\left( 1\right) x^{ss}=c \\
x^{ss} &=&\frac{c}{P\left( 1\right) }
\end{eqnarray}

provided $P\left( 1\right) \neq 0.$

In the more general case we often have to guess a particular solution.

\subsection{Systems of linear first order difference equations}

Systems of first order difference equations are solved with the
diagonalization method that we also used for the differential equation.

\begin{eqnarray}
\left[ 
\begin{array}{c}
x_{1,t+1} \\ 
x_{2,t+1} \\ 
. \\ 
x_{n,t+1}%
\end{array}
\right] &=&\mathbf{A}\left[ 
\begin{array}{c}
x_{1,t} \\ 
x_{2,t} \\ 
. \\ 
x_{n,t}%
\end{array}
\right] +\mathbf{P} \\
\mathbf{x}_{t+1} &=&\mathbf{Ax}_{t}+\mathbf{P.}
\end{eqnarray}

First we find the diagonalizing matrix of eigenvectors $\mathbf{B}^{-1}$.
Then we solve the homogeneous equation by defining 
\begin{eqnarray}
\mathbf{y}_{t+1} &=&\mathbf{Bx}_{t+1} \\
&=&\mathbf{BAx}_{t} \\
&=&\mathbf{BA\mathbf{B}^{-1}y}_{t} \\
&=&\left[ 
\begin{array}{cccc}
r_{1}^{t} & 0 & . & 0 \\ 
0 & r_{2}^{t} & . & 0 \\ 
. & . & . & . \\ 
0 & 0 & . & r_{n}^{t}%
\end{array}%
\right] \left[ 
\begin{array}{c}
c_{1} \\ 
c_{2} \\ 
. \\ 
c_{n}%
\end{array}%
\right] \\
&\equiv &\mathbf{r}^{t}\mathbf{c.}
\end{eqnarray}

Then we transform back\footnote{%
Since the vector of constants, $\mathbf{c}$, is arbitrary, we may equally
well write $\mathbf{x}_{t}=\mathbf{B}^{-1}\mathbf{r}^{t}\mathbf{c}$.} 
\begin{equation}
\mathbf{x}_{t+1}=\mathbf{B}^{-1}\mathbf{y}_{t+1}=\mathbf{B}^{-1}\mathbf{r}%
^{t}\mathbf{c.}  \label{Diff1}
\end{equation}

A particular solution to the complete equation then has to be added. Here we
might try to find a steady state as a particular solution

\begin{eqnarray}
\mathbf{x}^{ss} &=&\mathbf{Ax}^{ss}+\mathbf{P,} \\
\mathbf{x}^{ss} &=&\left( \mathbf{I-A}\right) ^{-1}\mathbf{P.}
\end{eqnarray}

If we have the initial conditions we find that 
\begin{eqnarray}
\mathbf{x}_{t+1} &=&\mathbf{B}^{-1}\mathbf{r}^{t}\mathbf{c+x}^{ss} \\
\mathbf{x}_{1} &=&\mathbf{B}^{-1}\mathbf{r}^{0}\mathbf{c+x}^{ss} \\
&=&\mathbf{B}^{-1}\mathbf{c+x}^{ss} \\
\mathbf{c} &=&\mathbf{B}\left( \mathbf{x}_{1}-\mathbf{x}^{ss}\right) \\
\mathbf{x}_{t+1} &=&\mathbf{B}^{-1}\mathbf{r}^{t}\mathbf{B}\left( \mathbf{x}%
_{1}-\mathbf{x}^{ss}\right) \mathbf{+x}^{ss}
\end{eqnarray}

An example;

\begin{eqnarray}
\mathbf{x}_{t+1} &=&\left[ 
\begin{array}{cc}
0 & 1 \\ 
-1 & 0%
\end{array}%
\right] \mathbf{x}_{t}+\left[ 
\begin{array}{c}
-1 \\ 
-1%
\end{array}%
\right] . \\
\mathbf{x}^{ss} &=&\left( \left[ 
\begin{array}{cc}
1 & 0 \\ 
0 & 1%
\end{array}%
\right] -\left[ 
\begin{array}{cc}
0 & 1 \\ 
-1 & 0%
\end{array}%
\right] \right) ^{-1}\left[ 
\begin{array}{c}
-1 \\ 
-1%
\end{array}%
\right] =\left[ 
\begin{array}{c}
-1 \\ 
0%
\end{array}%
\right]
\end{eqnarray}

Eigenvalues are 
\begin{eqnarray}
r_{1,2} &=&i,-i\rightarrow r_{1,2}^{t}=1^{t}\left( \cos \left( \frac{\pi }{2}%
t\right) \pm i\sin \left( \frac{\pi }{2}t\right) \right) \\
&=&\left( \cos \left( \frac{\pi }{2}t\right) \pm i\sin \left( \frac{\pi }{2}%
t\right) \right)
\end{eqnarray}
and 
\begin{equation}
\mathbf{B}^{-1}=\left[ 
\begin{array}{cc}
1 & 1 \\ 
i & -i%
\end{array}
\right] .
\end{equation}

The solution is 
\begin{gather}
\mathbf{x}_{t}=\left[ 
\begin{array}{cc}
1 & 1 \\ 
i & -i%
\end{array}%
\right] \left[ 
\begin{array}{cc}
i^{t} & 0 \\ 
0 & \left( -i\right) ^{t}%
\end{array}%
\right] \left[ 
\begin{array}{c}
c_{1} \\ 
c_{2}%
\end{array}%
\right] +\left[ 
\begin{array}{c}
-1 \\ 
0%
\end{array}%
\right] \\
=\left[ 
\begin{array}{cc}
i^{t} & \left( -i\right) ^{t} \\ 
i^{t+1} & \left( -i\right) ^{t+1}%
\end{array}%
\right] \left[ 
\begin{array}{c}
c_{1} \\ 
c_{2}%
\end{array}%
\right] +\left[ 
\begin{array}{c}
-1 \\ 
0%
\end{array}%
\right] \\
=\left[ 
\begin{array}{c}
c_{1}i^{t}+c_{2}\left( -i\right) ^{t} \\ 
c_{1}i^{t+1}+c_{2}\left( -i\right) ^{t+1}%
\end{array}%
\right] +\left[ 
\begin{array}{c}
-1 \\ 
0%
\end{array}%
\right] \\
=\left[ 
\begin{array}{c}
c_{1}\left( \left( \cos \left( \frac{\pi }{2}t\right) +i\sin \left( \frac{%
\pi }{2}t\right) \right) \right) \\ 
+c_{2}\left( \cos \left( \frac{\pi }{2}t\right) -i\sin \left( \frac{\pi }{2}%
t\right) \right) \\ 
c_{1}\left( \left( \cos \left( \frac{\pi }{2}\left( t+1\right) \right)
+i\sin \left( \frac{\pi }{2}\left( t+1\right) \right) \right) \right) \\ 
+c_{2}\left( \cos \left( \frac{\pi }{2}\left( t+1\right) \right) -i\sin
\left( \frac{\pi }{2}\left( t+1\right) \right) \right)%
\end{array}%
\right] \\
+\left[ 
\begin{array}{c}
-1 \\ 
0%
\end{array}%
\right] \\
=\left[ 
\begin{array}{c}
\left( c_{1}+c_{2}\right) \left( \cos \left( \frac{\pi }{2}t\right) +\left(
c_{1}-c_{2}\right) i\sin \left( \frac{\pi }{2}t\right) \right) \\ 
\left( c_{1}+c_{2}\right) \left( \cos \left( \frac{\pi }{2}\left( t+1\right)
\right) +\left( c_{1}-c_{2}\right) i\sin \left( \frac{\pi }{2}\left(
t+1\right) \right) \right)%
\end{array}%
\right] +\left[ 
\begin{array}{c}
-1 \\ 
0%
\end{array}%
\right] \\
\left[ 
\begin{array}{c}
\tilde{c}_{1}\left( \cos \left( \frac{\pi }{2}t\right) +\tilde{c}_{2}\sin
\left( \frac{\pi }{2}t\right) \right) \\ 
\tilde{c}_{1}\left( \cos \left( \frac{\pi }{2}\left( t+1\right) \right) +%
\tilde{c}_{2}\sin \left( \frac{\pi }{2}\left( t+1\right) \right) \right)%
\end{array}%
\right] +\left[ 
\begin{array}{c}
-1 \\ 
0%
\end{array}%
\right] .
\end{gather}

\subsubsection{Non-invertible eigenvectors}

If some roots are repeated $\mathbf{B}^{-1}$ may be non-invertible. In this
case we cannot use the diagonalization method. Instead we can use the
existence of a higher order single difference equation that is equivalent to
the system of first order difference equations we want to solve. This is
exactly analogous to the case of differential equations.

Look at the following example: 
\begin{equation}
\mathbf{x}_{t+1}=\left[ 
\begin{array}{cc}
-1 & 1 \\ 
-4 & 3%
\end{array}%
\right] \mathbf{x}_{t}  \label{differencesystem1}
\end{equation}

The eigenvalues of the coefficient matrix are both 1 and the matrix of
eigenvectors are non-invertible. Now we transform the system into a second
order single difference equation.

From the first row we have that

\begin{eqnarray}
x_{1,t+1} &=&-x_{1,t}+x_{2,t}  \label{subst} \\
&\rightarrow &x_{2,t}=x_{1,t+1}+x_{1,t}, \\
x_{2,t+1} &=&x_{1,t+2}+x_{1,t+1}
\end{eqnarray}

Using the second row, 
\begin{eqnarray}
x_{2,t+1} &=&-4x_{1,t}+3x_{2,t}  \label{equivsystem} \\
&\rightarrow &x_{1,t+2}+x_{1,t+1}=-4x_{1,t}+3\left( x_{1,t+1}+x_{1,t}\right)
\\
0 &=&x_{1,t+2}-2x_{1,t+1}+x_{1,t}
\end{eqnarray}

Note that the polynomial associated with this second order difference
equation is identical to the characteristic equation of the coefficient
matrix in (\ref{differencesystem1}). Consequently, they have the same
(repeated) root 1. The solution to (\ref{equivsystem}) is 
\begin{equation}
x_{1,t}=\left( c_{1}+tc_{2}\right) 1^{t}=\left( c_{1}+tc_{2}\right)
\end{equation}

With knowledge of $x_{1,0}$ and $x_{1,1}$, we have 
\begin{equation}
x_{1,t}=\left( x_{1,0}+t\left( x_{1,1}-x_{1,0}\right) \right)
\end{equation}

By using (\ref{subst}), we can express the solution in terms of $x_{1,t}$
and $x_{2,t}$%
\begin{eqnarray}
x_{1,t} &=&x_{1,0}+t\left( x_{1,1}-x_{1,0}\right) \\
x_{2,t} &=&x_{1,t+1}+x_{1,t} \\
&=&2x_{1,0}+\left( 2t+1\right) \left( x_{1,1}-x_{1,0}\right) \\
&\rightarrow &\left[ 
\begin{array}{c}
x_{1,t} \\ 
x_{2,t}%
\end{array}%
\right] =\left[ 
\begin{array}{cc}
1 & t \\ 
2 & 2t+1%
\end{array}%
\right] \left[ 
\begin{array}{c}
x_{1,0} \\ 
x_{1,1}-x_{1,0}%
\end{array}%
\right] .
\end{eqnarray}%
\newpage

\section{Dynamic Optimization in Discrete Time}

\subsection{Non-Stochastic Dynamic Programming}

%TCIMACRO{\TeXButton{TeX field}{\setcounter{equation}{0}}}%
%BeginExpansion
\setcounter{equation}{0}%
%EndExpansion
Consider the dynamic problem

\begin{eqnarray}
&&\max_{\left\{ k_{t},c_{t}\right\} _{t=1}^{T}}\sum_{t=1}^{T}u\left(
k_{t},c_{t},t\right)  \label{Dynoptproblem1} \\
\text{s.t.\thinspace }k_{1} &=&\underline{k}  \label{constraint1} \\
k_{t+1} &=&f\left( k_{t},c_{t},t\right) ,t=1...,T, \\
k_{T+1} &=&\bar{k}  \label{constraint3} \\
c_{t} &\in &\Omega _{t}
\end{eqnarray}

Before trying to solve this, notice

\begin{enumerate}
\item Per period payoff is additive over time.

\item $k_{t}$ cannot be changed in period $t$, but its future values, its
law-of-motion can be changed by $c_{t}.$ We will call $k$ a \emph{state
variable} (to be more properly defined later) and $c$ a \emph{control
variable }that must be chosen in the control set $\Omega _{t}$. A sequence $%
k_{1},k_{2},...,k_{T+1}$ is said to be \emph{admissible} if and only if it
satisfies the constraints (\ref{constraint1}-\ref{constraint3}) for some
sequence $\left\{ c_{t}\in \Omega _{t}\right\} _{t=1}^{T}$.
\end{enumerate}

A direct way to solve this would be to form the Lagrangian

\begin{eqnarray}
L &=&\sum_{t=1}^{T}u\left( k_{t},c_{t},t\right) +\sum_{t=1}^{T}\lambda
_{t}\left( f\left( k_{t},c_{t},t\right) -k_{t+1}\right) \\
&&+\mu _{1}\left( \underline{k}-k_{1}\right) +\mu _{T+1}\left( k_{T+1}-\bar{k%
}\right)
\end{eqnarray}

with first order conditions

\begin{eqnarray}
u_{k}\left( k_{t},c_{t},t\right) +\lambda _{t}f_{k}\left(
k_{t},c_{t},t\right) -\lambda _{t-1} &=&0,\forall t=2,..,T, \\
u_{c}\left( k_{t},c_{t},t\right) +\lambda _{t}f_{c}\left(
k_{t},c_{t},t\right) &=&0,\forall t=1,..,T, \\
u_{k}\left( k_{1},c_{1},1\right) +\lambda _{1}f_{k}\left(
k_{1},c_{1},1\right) -\mu _{1} &=&0, \\
-\lambda _{T}+\mu _{T+1} &=&0,
\end{eqnarray}

and (\ref{constraint1}-\ref{constraint3}).

This works, at least in principle, if $T$ is finite. An alternative way is
to recognize that in a problem like this, each sub-section of the path must
be optimal in itself. This means that the problem has a recursive
formulation i.e., it can be set up sequentially. We can thus solve the
problem backwards starting from the last period. In any period, the
remaining problem only depends on earlier actions through the ''inherited''
value of $k$.

For example, if the problem is over three periods ($T=3$) we can rewrite (%
\ref{Dynoptproblem1}) 
\begin{gather}
\max_{c_{1},k_{2}|_{k_{1}}}\left( u\left( k_{1},c_{1},1\right)
+\max_{c_{2},k_{3}|_{k_{2}}}\left( u\left( k_{2},c_{2},2\right)
+\max_{c_{3},k_{4}|_{k_{3}}}u\left( k_{3},c_{3},3\right) \right) \right) \\
\text{s.t. }k_{t+1}=f\left( k_{t},c_{t},t\right) ,t=1...,3 \\
k_{4}=\bar{k}
\end{gather}

In the final period ($T=3$), the problem is trivial; simply set $c_{3}$ so
that $k_{4}=\bar{k}$. The value of $c_{3}$ that solves $\bar{k}=f\left(
k_{3},c_{3},3\right) $ is a function of $k_{3}$.\footnote{%
For now, we just assume it is unique.} Denote that function $c_{3}\left(
k_{3}\right) .$ We can then define 
\begin{equation}
u\left( k_{3},c_{3}\left( k_{3}\right) ,3\right) \equiv W\left(
k_{3},3\right) .
\end{equation}

The interpretation of $W\left( k_{3},3\right) $, is the \emph{maximum
remaining pay-off in period 3,} being a function of the state variable $%
k_{3} $.

In period 2, we then want to solve 
\begin{eqnarray}
&&\max_{c_{2},k_{3}}\left( u\left( k_{2},c_{2},2\right) +u\left(
k_{3},c_{3}\left( k_{3}\right) ,3\right) \right) \\
\text{s.t.\thinspace }k_{3} &=&f\left( k_{2},c_{2},2\right) .
\end{eqnarray}

Using $W\left( k_{3},3\right) ,$we can write this 
\begin{equation}
\max_{c_{2}}\left( u\left( k_{2},c_{2},2\right) +W\left( f\left(
k_{2},c_{2},2\right) ,3\right) \right)
\end{equation}%
The solution and the maximized value depends on $k_{2}$ only and we called
the latter the value function and $k_{2}$ the state variable. We can then
define 
\begin{equation}
W\left( k_{2},2\right) \equiv \max_{c_{2}}\left( u\left(
k_{2},c_{2},2\right) +W\left( f\left( k_{2},c_{2},2\right) ,3\right) \right)
.
\end{equation}

The interpretation of this Bellman equation is straightforward. It says that
the maximum remaining pay-off in period 2, being a function of $k_{2}$, is
identically (i.e., for all $k_{2}$) equal to the maximum over the control in
period 2, $c_{2},$ over period 2 pay-off and the maximum remaining pay-off
in period 3 with period 3 state variable given by $f\left(
k_{2},c_{2},2\right) .$

The trade-off between generation of current and future pay-off is optimized
by \emph{one} simple FOC 
\begin{equation}
u_{c}\left( k_{2},c_{2},2\right) +W_{k}\left( f\left( k_{2},c_{2},2\right)
,3\right) f_{c}\left( k_{2},c_{2},2\right) =0.
\end{equation}

Finally, in the first period, 
\begin{equation}
W\left( k_{1},1\right) \equiv \max_{c_{1}}\left( u\left(
k_{1},c_{1},1\right) +W\left( f\left( k_{1},c_{1},1\right) ,2\right) \right)
.
\end{equation}

If we know the value functions, the multidimensional problem has become much
simpler. The Bellman equation provides a way of verifying that the value
function we use is correct. It is of course straightforward to extend the
analysis to any finite horizon problem, yielding 
\begin{eqnarray}
W\left( k_{t},t\right) &\equiv &\max_{c_{t}}\left( u\left(
k_{t},c_{t},t\right) +W\left( k_{t+1},t+1\right) \right) \\
\text{s.t. }k_{t+1} &=&f\left( k_{t},c_{t},t\right)
\end{eqnarray}

\subsubsection{Discounting and the Current Value Bellman equation}

Very often in macroeconomics, the objective function is a discounted sum of
pay-offs, i.e., (\ref{Dynoptproblem1}) can be written%
\begin{equation}
\max_{\left\{ c\right\} _{t=1}^{T}}\sum_{t=1}^{T}\beta ^{t-1}u\left(
k_{t},c_{t},t\right) .
\end{equation}%
In this case, it is convenient to work with \emph{current value functions, }$%
V\left( k,t\right) $ which should be interpreted as the maximum remaining
value that can be achived from time $t$ and onward, given $k_{t}$ and \emph{%
seen from period }$t.$ In other words, given a problem 
\begin{eqnarray}
&&\max_{\left\{ c\right\} _{t=1}^{T}}\sum_{t=1}^{T}\beta ^{t-1}u\left(
k_{t},c_{t},t\right) \\
\text{s.t.~}k_{t+1} &=&f\left( k_{t},c_{t},t\right) ,t=1...,T \\
k_{1} &=&\underline{k} \\
k_{T+1} &=&\bar{k}  \notag
\end{eqnarray}%
we define for any $t\in \left\{ 1,...,T\right\} $%
\begin{eqnarray}
V\left( k_{t},t\right) &\equiv &\max_{\left\{ c,k\right\}
_{s=t}^{T}}\sum_{s=t}^{T}\beta ^{s-t}u\left( k_{s},c_{s},s\right) \\
\text{s.t.~}k_{t+1} &=&f\left( k_{t},c_{t},t\right) ,t=s...,T. \\
k_{T+1} &=&\bar{k}.  \notag
\end{eqnarray}

From this follows that 
\begin{equation*}
V\left( k,t\right) \equiv \beta ^{-\left( t-1\right) }W\left( k,t\right) .
\end{equation*}

Using this in the Bellman equation, (where I substitute for $k_{t+1}$ from
the law-of-motion) we get the \emph{current value Bellman equation}%
\begin{eqnarray}
W\left( k_{t},t\right) &\equiv &\max_{c_{t}}\left( \beta ^{t-1}u\left(
k_{t},c_{t},t\right) +W\left( f\left( k_{t},c_{t},t\right) ,t+1\right)
\right) \\
\beta ^{t-1}V\left( k,t\right) &\equiv &\max_{c_{t}}\left( \beta
^{t-1}u\left( k_{t},c_{t},t\right) +\beta ^{t}V\left( f\left(
k_{t},c_{t},t\right) ,t+1\right) \right) \\
V\left( k,t\right) &\equiv &\max_{c_{t}}\left( u\left( k_{t},c_{t},t\right)
+\beta V\left( f\left( k_{t},c_{t},t\right) ,t+1\right) \right)
\label{CurrentBellman}
\end{eqnarray}%
In practice, the current value Bellman equation is the most used variant in
macroeconomics and, therefore, you will often see the word \emph{current }%
dropped \ and (\ref{CurrentBellman}) is simple referred to as the Bellman
equation and $V\left( k,t\right) $ is referred to as the value function.

\subsubsection{Infinite Horizon and Autonomous Problems}

In an infinite horizon problem we cannot use the method of starting from the
last period. Still, if the problem has a well-defined value function, \emph{%
it satisfies the Bellman equation}. Furthermore, under some conditions
(which we will return to later), there is\emph{\ only one} function that
solves the Bellman equation, so if we find one function that solves the
Bellman equation, we have a solution to the dynamic optimization problem.
Since geometric discounting will prove to be important for showing
uniqueness, we will use that from now on. An important case when there is
necessarily exists a unique solution to Bellman equation is if there is
strict discounting ($\beta \in \left( 0,1\right) $) and $u\left(
k_{t},c_{t},t\right) $ is bounded for any admissible sequence of $\left\{
k_{t},c_{t}\right\} .$

To find a solution, we will use two different approaches.

\begin{enumerate}
\item Guess on a value function and verify that it satisfies the Bellman
equation.

\item Iterate on the Bellman equation until it converges.
\end{enumerate}

Guessing is often feasible when the problem is autonomous (stationary).
Then, the problem is independent of time in the sense that given an initial
condition on the state variable(s), the solution and the maximized objective
is independent of the starting date. A problem is autonomous if

\begin{enumerate}
\item Time is infinite,

\item the law of motion for the state (including constraints on the control)
is independent of time, and

\item the per-period return function is the same over time, except for
possibly a geometric discount factor, i.e., $u\left( k,c,t\right) =\beta
^{t}u\left( k,c\right) $.
\end{enumerate}

(Think about what would happen if any of these conditions is not satisfied).
In this case, the \emph{current }value function turns out to be independent
of time.\footnote{%
The present value function $W(k,t)$ is not independent ot time, but
separable so that we can write $W\left( k,t\right) =W\left( k\right) \beta
^{t-1}.$} We can then write the current Bellman equation in terms of the
current value function 
\begin{eqnarray}
V\left( k_{t}\right) &\equiv &\max_{c_{t}}\left( u\left( k_{t},c_{t}\right)
+\beta V\left( k_{t+1}\right) \right)  \label{Bellman1} \\
\text{s.t. }k_{t+1} &=&f\left( k_{t},c_{t}\right) .
\end{eqnarray}

Suppose we find a solution to the maximization problem in the RHS\ of (\ref%
{Bellman1}). This will be a time invariant function $c(k)$, since $u,$ $V$
and $f$ are time-invariant. Plugging $c\left( k\right) $ into (\ref{Bellman1}%
), we get rid of the max-operator: 
\begin{equation}
V\left( k\right) =u\left( k,c\left( k\right) \right) +\beta V\left( f\left(
k,c\left( k\right) \right) \right)  \label{Bellmanidentity1}
\end{equation}

If (\ref{Bellmanidentity1}) is satisfied for \emph{all values of admissible }%
$\emph{k,}$ we have a solution to the value function, otherwise our guess
was incorrect.

Note that (\ref{Bellmanidentity1}) is a functional equation, i.e., the LHS\
and RHS have to be the same functions. It is convenient to define the RHS\
as 
\begin{equation}
T(V\left( k\right) )\equiv \max_{c}u\left( k,c\right) +\beta V\left( f\left(
k,c\right) \right)
\end{equation}%
where $T$ operates on \emph{functions} rather than on values.\footnote{%
See section \ref{Tmap} for more on the $T$-mapping.} In the autonomous case,
when the value function is unchanged over time, the Bellman equation defines
a \emph{fixed point} for $T$ in the space of functions $V$; 
\begin{equation}
V\left( k\right) =T\left( V\left( k\right) \right) .  \label{Fixedpoint}
\end{equation}

(\ref{Fixedpoint}) means that if we plug in some function of $k$ in the RHS
of (\ref{Fixedpoint}) we must get out the same function on the LHS. The
Bellman equation is a necessary condition for $V\left( k\right) $ being a
correctly specified value function, we will later discuss conditions under
which it is also sufficient.

Typically the value function and the objective function are of similar form.
This is intuitive in the light of (\ref{Bellmanidentity1}). For example, if
the $u$ function is logarithmic, we guess that the value function is of the
logarithmic form too. For HARA utility functions (e.g., CRRA, CARA and
quadratic) the value functions are generally of the same type as the utility
function (Merton, 1971).

\subsubsection{An example of guessing}

In the problem (\ref{Dynoptproblem1}) let time be infinite and 
\begin{eqnarray}
u(k,c,t) &=&\beta ^{t}\ln (c), \\
f(k,c,t) &=&k^{\alpha }-c,0<\alpha <1,
\end{eqnarray}
and let the end-condition $k_{T+1}=\bar{k}$ be replaced by $k_{t}>0$ $%
\forall t.$

This is an autonomous problem, so we have 
\begin{eqnarray}
V\left( k_{t}\right) &=&\max_{c,k}\left( u\left( c_{t}\right) +\beta V\left(
k_{t+1}\right) \right) \\
\text{s.t. }k_{t+1} &=&f\left( k_{t},c_{t}\right) , \\
&\Rightarrow &V\left( k_{t}\right) =\max_{c}\left( \ln c_{t}+\beta V\left(
k_{t}^{\alpha }-c_{t}\right) \right)
\end{eqnarray}

Now, guess that $V$ is of the same form as $u$, here $X\ln k_{t}+Y$, for
some unknown constants $X$ and $Y,$ giving first order conditions 
\begin{eqnarray}
u^{\prime }\left( c_{t}\right) &=&\beta V^{\prime }\left( k_{t}^{\alpha
}-c_{t}\right) \\
\frac{1}{c_{t}} &=&\beta \frac{X}{k_{t}^{\alpha }-c_{t}} \\
&\Rightarrow &c_{t}=\frac{k_{t}^{\alpha }}{1+\beta X}\equiv c\left(
k_{t}\right) . \\
k_{t+1} &=&k_{t}^{\alpha }-c_{t}=k_{t}^{\alpha }-\frac{k_{t}^{\alpha }}{%
1+\beta X}=\frac{\beta X}{1+\beta X}k_{t}^{\alpha }
\end{eqnarray}

Plugging $c\left( k_{t}\right) $ into the Bellman equation yields 
\begin{eqnarray}
X\ln k_{t}+Y &=&\ln \frac{k_{t}^{\alpha }}{1+\beta X}+\beta \left( X\ln 
\frac{\beta X}{1+\beta X}k_{t}^{\alpha }+Y\right) \\
&=&\ln \frac{k_{t}^{\alpha }}{1+\beta X}+\beta \left( X\ln \frac{%
k_{t}^{\alpha }}{1+\beta X}+X\ln \beta X+Y\right) \\
&=&\left( 1+\beta X\right) \ln \frac{k_{t}^{\alpha }}{1+\beta X}+\beta
\left( X\ln \beta X+Y\right) \\
&=&\alpha \left( 1+\beta X\right) \ln k_{t}-\left( 1+\beta X\right) \ln
\left( 1+\beta X\right) \\
&&+\beta X\ln \beta X+\beta Y.  \notag
\end{eqnarray}

This is true for all values of $k$, if and only if 
\begin{eqnarray}
X &=&\alpha \left( 1+\beta X\right) \\
Y &=&-\left( 1+\beta X\right) \ln \left( 1+\beta X\right) +\beta X\ln \beta
X+\beta Y.
\end{eqnarray}
giving 
\begin{eqnarray}
X &=&\frac{\alpha }{1-\alpha \beta } \\
Y &=&\frac{\beta \alpha \ln \left( \beta \alpha \right) +\left( 1-\beta
\alpha \right) \ln \left( 1-\alpha \beta \right) }{\left( 1-\beta \right)
\left( 1-\alpha \beta \right) } \\
&\Rightarrow &V\left( k\right) =\frac{\alpha }{1-\alpha \beta }\ln k+\frac{%
\beta \alpha \ln \left( \beta \alpha \right) +\left( 1-\beta \alpha \right)
\ln \left( 1-\alpha \beta \right) }{\left( 1-\beta \right) \left( 1-\alpha
\beta \right) }.
\end{eqnarray}

Having $V\left( k\right) ,$ it is easy to find the \emph{optimal control},
or the \emph{policy rule}, 
\begin{eqnarray}
c\left( k_{t}\right) &=&\frac{k_{t}^{\alpha }}{1+\beta \frac{\alpha }{%
1-\alpha \beta }}=\left( 1-\alpha \beta \right) k_{t}^{\alpha }, \\
\ln k_{t+1} &=&\ln \alpha \beta +\alpha \ln k_{t.}
\end{eqnarray}

\subsubsection{Iteration}

As mentioned above, an alternative way to find the infinite horizon value
function is to find the limit of finite horizon Bellman equation as the
horizon goes to infinity. Under for our purposes quite general conditions
this limit exists and is equal to the value function for the infinite
horizon problem. Let us change notation slightly, measuring time as the
number of remaining periods $s$ until the final period. We then denote the
(current) value function with $s$ periods left by 
\begin{equation}
V(k,s).
\end{equation}%
\qquad\ 

We also assume geometric discounting and that both pay-offs and the
law-of-motion for the state variable are time-independent ($u(k,c,s)=u\left(
k,c\right) $ and $f\left( k,c,s\right) =f\left( k,c\right) $) so that the
infinite horizon problem is autonomous. If the following limit is
well-defined, we denote 
\begin{equation}
\lim {}_{s\rightarrow \infty }V(k,s)\equiv V\left( k\right) .
\end{equation}

The iteration method is usually done numerically, but it can (at some cost
of messiness) be done also analytically. Using the $T$ operator and the
Bellman equation, we find the limit in the following way 
\begin{eqnarray}
V\left( k,s\right) &=&\max_{c_{s}}\left( u\left( k_{s},c_{s}\right) +\beta
V\left( f\left( k_{s},c_{s}\right) ,s-1\right) \right) \\
&\equiv &T\left( V\left( k,s-1\right) \right) \\
V\left( k,s\right) &=&T^{s}V\left( k,0\right) \\
V\left( k\right) &\equiv &\lim_{s\rightarrow \infty }V\left( k,s\right)
=\lim_{s\rightarrow \infty }T^{s}V\left( k,0\right) .
\end{eqnarray}

If the limit exists, it clearly satisfies the Bellman equation 
\begin{eqnarray}
V\left( k\right) &=&T\left( V\left( k\right) \right) \\
\lim_{s\rightarrow \infty }T^{s}V\left( k,0\right) &=&T\left(
\lim_{s\rightarrow \infty }T^{s}V\left( k,0\right) \right)
=\lim_{s\rightarrow \infty }T^{s+1}V\left( k,0\right) .
\end{eqnarray}

The remaining issue is what function $V\left( k,0\right) $ to start the
iteration with. However, suppose that we can show that the limit $V\left(
k\right) $ satisfies 
\begin{equation}
\lim_{s\rightarrow \infty }\beta ^{s}V\left( k_{s}\right) =0
\label{iterationlimit}
\end{equation}%
for ALL\ admissable\ values of $k_{t+s}$ that can be reached, given relevant
initial conditions and other constraints. Then, it can easily be shown that
the Bellman equation is \emph{sufficient. }Then, $\lim_{s\rightarrow \infty
}T^{s}V\left( k,0\right) $ provides a unique solution to the Bellman
equation. This means that the limit is independent of the choice of $V\left(
k,0\right) .$ As we see from (\ref{iterationlimit}), $\beta <1$ and $V\left(
k\right) $ bounded are sufficient for uniqueness. Intuitively, if $\beta <1$
and pay-offs are bounded, the pay-off in the infinite horizon has no impact
on the value function. Let us revert to measure time in the usual way. Then,
if the Bellman equation is satisfied, we have 
\begin{gather}
V\left( k_{t}\right) =\max_{c_{t}}\left( u\left( k_{t},c_{t}\right) +\beta
V\left( k_{t+1}\right) \right) \\
\text{s.t. }k_{t+1}=f\left( k_{t},c_{t}\right)  \notag \\
\max_{c_{t}}\left( u\left( k_{t},c_{t}\right) +\beta \left(
\max_{c_{t+1}}\left( u\left( k_{t+1},c_{t+1}\right) +\beta V\left(
k_{t+2}\right) \right) \right) \right) \\
\text{s.t.}k_{t+1}=f\left( k_{t},c_{t}\right) ,k_{t+2}=f\left(
k_{t+1},c_{t+1}\right)  \notag \\
=\max_{\{c_{t+n}\}_{0}^{1}}\sum_{n=0}^{1}\beta ^{n}u\left(
k_{t+n},c_{t+n}\right) +\beta ^{2}V\left( f\left( k_{t+1},c_{t+1}\right)
\right) , \\
\text{s.t.}k_{t+1}=f\left( k_{t},c_{t}\right) ,k_{t+2}=f\left(
k_{t+1},c_{t+1}\right)
\end{gather}

Repeating this, and taking the limit yields 
\begin{eqnarray}
V\left( k_{t}\right) &=&\max_{\{c_{t+n}\}_{0}^{s}}\sum_{n=0}^{s}\beta
^{n}u\left( k_{t+n},c_{t+n}\right) +\beta ^{s+1}V\left( k_{t+n+1}\right) , \\
\text{s.t.}k_{t+n+1} &=&f\left( k_{t+n},c_{t+n}\right) \forall n\in \left\{
0,s\right\} \\
V\left( k_{t}\right) &=&\max_{\{c_{t+n}\}_{0}^{\infty }}\sum_{n=0}^{\infty
}\beta ^{n}u\left( k_{t+n},c_{t+n}\right) +\lim_{s\rightarrow \infty }\beta
^{s+1}V\left( k_{t+s+1}\right) , \\
\text{s.t.}k_{t+n+1} &=&f\left( k_{t+n},c_{t+n}\right) \forall n\geq 0,
\end{eqnarray}%
Now, if $\lim_{s\rightarrow \infty }\beta ^{s+1}V\left( k_{t+s+1}\right) =0,$
for all permissible paths of $k$ we have showed that 
\begin{eqnarray}
V\left( k_{t}\right) &=&\max_{\{c_{t+s}\}_{0}^{\infty }}\sum_{s=0}^{\infty
}\beta ^{s}u\left( k_{t+s},c_{t+s}\right) \\
\text{s.t. }k_{t+s+1} &=&f\left( k_{t+s},c_{t+s}\right) \forall s\geq 0\text{%
, given}k_{t},
\end{eqnarray}%
i.e., that the Bellman equation implies optimality.

The iteration can easily be done numerically, either by specifying a
functional form if we know that, or by just choosing a grid. In the latter
case we assume that the state variable must belong to a finite set of
values, say for example that in every period $k$ must be chosen from the set 
$\mathbf{K\equiv }\left\{ k_{1},k_{2},...,k_{n}\right\} $. Then, we can
compute the corresponding set of possible controls, $c_{m,n}\in \mathbf{C}$
from the equation 
\begin{equation}
k_{m}=f\left( k_{n},c_{m,n}\right) .
\end{equation}

Then, in each iteration, we solve the Bellman equation for each $k\in 
\mathbf{K}$, giving for iteration $s$%
\begin{equation}
V\left( k_{n},s\right) =\max_{c_{m,n}\in \mathbf{C}}\left( u\left(
k_{n},c_{m,n}\right) +\beta V\left( k_{m},s-1\right) \right) .
\end{equation}

This goes quickly on a computer and the iteration is repeated until $V\left(
k,s\right) $ is sufficiently close to $V\left( k,s-1\right) $ over the set
of $k\in \mathbf{K.}$

\subsubsection{An envelope result}

We will later have use for the following envelope result, which implies that
we evaluate $dV\left( k\right) /dk\,$as the partial derivative holding $c$
constant. To see this, note that 
\begin{eqnarray}
\frac{dV\left( k_{t}\right) }{dk_{t}} &\equiv &V^{\prime }\left(
k_{t}\right) =\frac{\partial u\left( k_{t},c_{t}\right) }{\partial k_{t}}%
+\beta V^{\prime }\left( k_{t+1}\right) \frac{\partial f\left(
k_{t},c_{t}\right) }{\partial k_{t}} \\
&&+\frac{dc_{t}}{dk_{t}}\left( \frac{\partial u\left( k_{t},c_{t}\right) }{%
\partial c_{t}}+\beta V^{\prime }\left( k_{t+1}\right) \frac{\partial
f\left( k_{t},c_{t}\right) }{\partial c_{t}}\right) .
\end{eqnarray}%
However, the Bellman equation implies that the last term above is zero if $c$
is chosen optimally, either because the first-order condition for an
interior maximum is satisfied $\frac{\partial u\left( k,c\right) }{\partial c%
}+\beta V^{\prime }\left( k_{t+1}\right) \frac{\partial f\left(
k_{t},c_{t}\right) }{\partial c_{t}}=0$, or in the case of a corner because $%
\frac{dc_{t}}{dk_{t}}=0.$

This envelope result is often very useful. Consider the example 
\begin{eqnarray}
u\left( k_{t},c_{t}\right) &=&v\left( c_{t}\right) , \\
f\left( k_{t},c_{t}\right) &=&y\left( k_{t}\right) -c_{t}.
\end{eqnarray}

The interior solution to the Bellman equation satisfies 
\begin{eqnarray}
0 &=&v^{\prime }\left( c_{t}\right) +\beta V^{\prime }\left( k_{t+1}\right) 
\frac{\partial f\left( k_{t},c_{t}\right) }{\partial c_{t}},  \label{consopt}
\\
&\rightarrow &v^{\prime }\left( c_{t}\right) =\beta V^{\prime }\left(
k_{t+1}\right) .
\end{eqnarray}
The envelope condition yields 
\begin{eqnarray}
V^{\prime }\left( k_{t}\right) &=&\frac{\partial u\left( k_{t},c_{t}\right) 
}{\partial k_{t}}+\beta V^{\prime }\left( k_{t+1}\right) \frac{\partial
f\left( k_{t},c_{t}\right) }{\partial k_{t}} \\
&=&v^{\prime }\left( c_{t}\right) y^{\prime }\left( k_{t}\right) . \\
V^{\prime }\left( k_{t+1}\right) &\rightarrow &v^{\prime }\left(
c_{t+1}\right) y^{\prime }\left( k_{t+1}\right) .
\end{eqnarray}
Using this in (\ref{consopt}) yields the \emph{Euler equation} 
\begin{equation}
v^{\prime }\left( c_{t}\right) =\beta v^{\prime }\left( c_{t+1}\right)
y^{\prime }\left( k_{t+1}\right) .
\end{equation}

\subsubsection{State Variables}

We often solve the dynamic programming problem by guessing a form of the
value function. The first thing to determine is then which variables should
enter, i.e., which variables are the state variables. The state variables
must satisfy both following conditions:

1. To enter the value function at time they must be realized at $t$.

Note, however, that it sometimes may be convenient to use a conditional
expectation $E_{t}(z_{t+s})$ as a state variable. The expectation as of $t$
is certainly realized at $t$ even if the stochastic variable $z_{t+s}$ is
not realized.

2. The set of variables chosen as state variables must together give
sufficient information so that the value of the program from $t$ and onwards
when the optimal control is chosen can be calculated.\footnote{%
Can you figure out what do we need if the per period utility function in (%
\ref{Dynoptproblem1}) were $u\left( c_{t},c_{t-1}\right) ?$}

Note, we should try to find the smallest such set. If, for example we have
an investment problem with several assets to invest in and without any costs
of adjusting the portfolio, total wealth may be a sufficient as a state
variable.

\subsection{Stochastic Dynamic Programming}

As long as the recursive structure of the problem is intact adding a
stochastic element to the transition equation does not change the Bellman
equation. Consider the problem 
\begin{eqnarray}
&&\max_{\{c_{t}\}_{0}^{\infty }}E_{0}\sum_{t=0}^{\infty }\beta ^{t}u\left(
k_{t},c_{t}\right) \\
\text{s.t. }k_{t+1} &=&f\left( k_{t},c_{t},\varepsilon _{t+1}\right) \forall
t\geq 0,k_{0}\,\text{given, and} \\
k_{t} &>&0\text{ }\forall t\text{.}
\end{eqnarray}%
where $E_{0}$ is the expectations operator, conditional on time $0$
information and we assume that $c_{t}$ can be chosen conditional on
information about $\varepsilon _{s}$ for all $s\leq t.$ Furthermore, let us
assume that the distribution of $\varepsilon _{t}$ is $i.i.d.$ over time.
Then, the Bellman equation becomes 
\begin{equation}
V\left( k_{t}\right) =\max_{c_{t}}\left( u\left( k_{t},c_{t}\right) +\beta
E_{t}V\left( f\left( k_{t},c_{t},\varepsilon _{t+1}\right) \right) \right) ,
\end{equation}%
with a first-order condition 
\begin{equation}
0=u_{c}\left( k_{t},c_{t}\right) +\beta E_{t}\left( V^{\prime }\left(
k_{t+1}\right) f_{c}\left( k_{t},c_{t},\varepsilon _{t+1}\right) \right) .
\end{equation}

Note that, in general $E_{t}\left( V^{\prime }\left( k_{t+1}\right)
f_{c}\left( k_{t},c_{t},\varepsilon _{t+1}\right) \right) \neq
E_{t}V^{\prime }\left( k_{t+1}\right) E_{t}f_{c}\left(
k_{t},c_{t},\varepsilon _{t+1}\right) $.

\subsubsection{A Stochastic Consumption Example}

Consider the following problem 
\begin{eqnarray}
&&\max_{\{c_{t}\}_{0}^{\infty }}E_{0}\sum_{t=0}^{\infty }\beta ^{t}\ln c_{t}
\\
\text{s.t. }A_{t+1} &=&\left( A_{t}-c_{t}\right) \left( 1+\tilde{m}%
_{t+1}\right) \forall t\geq 0, \\
A_{t} &\geq &0\text{, }\forall t\geq 0\text{, }A_{0}\,\text{given.}
\end{eqnarray}

The consumer decides how much to consume each period. The savings is placed
in a risky asset with gross return $\left( 1+\tilde{m}_{t+1}\right) $, that
is drawn from an i.i.d. distribution with $E\left( \ln \left( 1+\tilde{m}%
_{t+1}\right) \right) =$ $m>-\infty $. If, for example the gross return is
log-normal, with mean $\bar{m}$ and variance $\sigma ^{2}$ then, $E\left( 1+%
\tilde{m}_{t+1}\right) =e^{\bar{m}+\frac{\sigma ^{2}}{2}}.$

The problem is autonomous so we write the current value Bellman equation
with time independent value function $V$%
\begin{equation}
V\left( A_{t}\right) =\max_{c_{t}}\left\{ \ln c_{t}+\beta E_{t}V\left(
\left( A_{t}-c_{t}\right) \left( 1+\tilde{m}_{t+1}\right) \right) \right\} .
\end{equation}

The necessary first order condition for $c_{t}$ yield 
\begin{equation}
\frac{1}{c_{t}}=\beta E_{t}V^{^{\prime }}\left( A_{t+1}\right) \left( 1+%
\tilde{m}_{t+1}\right) .  \label{FOCStoch}
\end{equation}

Now we use Merton's result and guess that the value function is 
\begin{equation}
V\left( A_{t}\right) =Y+X\ln A_{t},
\end{equation}

for some constants $Y$ and $X$. Substituting into (\ref{FOCStoch}), we get 
\begin{eqnarray}
\frac{1}{c_{t}} &=&\beta E_{t}\frac{X}{A_{t+1}}\left( 1+\tilde{m}%
_{t+1}\right) , \\
&=&\beta E_{t}\frac{X\left( 1+\tilde{m}_{t+1}\right) }{\left(
A_{t}-c_{t}\right) \left( 1+\tilde{m}_{t+1}\right) }, \\
&=&\beta \frac{X}{A_{t}-c_{t}}, \\
&\rightarrow &c_{t}=\frac{A_{t}}{1+\beta X}, \\
A_{t}-c_{t} &=&\frac{\beta XA_{t}}{1+\beta X}.
\end{eqnarray}

Now we have to solve for the constant $X$. This is done by substituting the
solutions to the first order conditions and the guess into the Bellman
equation, 
\begin{gather}
Y+X\ln A_{t}=\max_{c_{t}}\left\{ \ln c_{t}+\beta E_{t}V\left( \left(
A_{t}-c_{t}\right) \left( 1+\tilde{m}_{t+1}\right) \right) \right\} \\
=\ln \frac{A_{t}}{1+\beta X}+\beta E_{t}V\left( \frac{\beta XA_{t}}{1+\beta X%
}\left( 1+\tilde{m}_{t+1}\right) \right) \\
=\ln \frac{A_{t}}{1+\beta X}+\beta E_{t}\left( Y+X\ln \left( \frac{\beta
XA_{t}}{1+\beta X}\left( 1+\tilde{m}_{t+1}\right) \right) \right) \\
=\left( 1+\beta X\right) \ln A_{t}-\left( 1+\beta X\right) \ln \left(
1+\beta X\right) +\beta Y \\
+\beta X\ln \beta X+\beta Xm.
\end{gather}

This is satisfied for all $A_{t}$ iff 
\begin{eqnarray}
X &=&\left( 1+\beta X\right) \\
&\rightarrow &X=\frac{1}{1-\beta }, \\
Y &=&-\frac{1}{1-\beta }\ln \left( \frac{1}{1-\beta }\right) +\beta Y+\beta 
\frac{1}{1-\beta }\ln \beta \frac{1}{1-\beta }+\beta \frac{m}{1-\beta } \\
&\rightarrow &Y=-\frac{1}{1-\beta }\ln \frac{1}{1-\beta }+\frac{\beta }{%
\left( 1-\beta \right) ^{2}}\left( m+\ln \beta \right) .
\end{eqnarray}

Thus, 
\begin{eqnarray}
c_{t} &=&\frac{A_{t}}{1+\beta \frac{1}{1-\beta }}=\left( 1-\beta \right)
A_{t}, \\
A_{t}-c_{t} &=&\beta A_{t}. \\
A_{t+1} &=&A_{t}\beta \left( 1+\tilde{m}_{t+1}\right) , \\
\ln A_{t+1} &=&\ln A_{t}+\ln \beta +\ln \left( 1+\tilde{m}_{t+1}\right) \\
E_{t}\ln A_{t+1} &=&\ln A_{t}+\ln \beta +m.
\end{eqnarray}%
Note that since $\ln \left( 1+\tilde{m}_{t+1}\right) $ is normally
distributed, $\left( 1+\tilde{m}_{t+1}\right) >0,$ for all $t,$ implying $%
A_{t}>0$ for all $t.$ If, on the other hand, $\left( 1+\tilde{m}%
_{t+1}\right) $ can be negative with positive probability, $E_{t}\ln \left(
1+\tilde{m}_{t+1}\right) $ is minus infinity implying that the value
function is ill-defined.

\subsection{Contraction mappings}

In the previous section we discussed guessing on solutions to the Bellman
equation. However, we would like to know whether there exists a solution and
whether it is unique. If the latter is not the case, it is not in principle
sufficient to guess and verify, since we might have other value functions
that also satisfy the Bellman equation. To prove existence and uniqueness we
will apply a contraction mapping argument.\footnote{%
An alternative is sometimes to look for the limit $\lim_{s\rightarrow \infty
}T^{s}\left( V(k,0)\right) $, which typically is the solution we are
interested in (at least in macroeconomics).} For this purpose, we first have
to define some concepts. Before, let's look at a simple example of a
contaction mapping. Below, I show a graph of Sweden. Some of you might
recognize the dots as representing some cities and towns in our country.

\FRAME{itbpFU}{3.1609in}{4.0049in}{0in}{\Qcb{Map of Sweden}}{}{slide13.emf}{%
\special{language "Scientific Word";type "GRAPHIC";maintain-aspect-ratio
TRUE;display "USEDEF";valid_file "F";width 3.1609in;height 4.0049in;depth
0in;original-width 10.0016in;original-height 7.5135in;cropleft
"0.2265";croptop "0.9495";cropright "0.6979";cropbottom "0.1510";filename
'Figures/Slide13.EMF';file-properties "XNPEU";}}

Next, let's put the same map in a smaller scale on top of the first. This is
done in the next picture.

\FRAME{itbpF}{3.1462in}{4.0058in}{0in}{}{}{slide14.emf}{\special{language
"Scientific Word";type "GRAPHIC";maintain-aspect-ratio TRUE;display
"USEDEF";valid_file "F";width 3.1462in;height 4.0058in;depth
0in;original-width 10.0016in;original-height 7.5135in;cropleft
"0.2270";croptop "0.9496";cropright "0.6972";cropbottom "0.1511";filename
'Figures/Slide14.EMF';file-properties "XNPEU";}}

If you look closely, you will see that two dots (cities) overlap. In fact,
they represent the nice city Karlstad in V\"{a}rmland. Perhaps more
importantly (at least for some), you cannot find any other dots overlapping.
In fact, you might be able to convince yourself that no other points in
Sweden than Karlstad overlap.

If you would move around the map but keeping the smaller map inside the
borders of the larger, you will be able to verify for yourself that there is
always one and one only points in Sweden that overlap in the two maps. One
can actually prove this result -- and generalizing this proof gives us the
contraction mapping proof.

Here is a sketch of how the result is formalized. Consider the space of
points in the large map. Put the smaller scale map ontop and this of this as
representing a mapping from the all points in large map back to the same set
of points. Karlstad maps back to Karlstad, Malm\"{o} maps to a point a
little west of J\"{o}nk\"{o}ping, Kiruna to a point northwest of Sundsvall
and so on.

Lets call this mapping $T$ and the set of points in the large map $X.$
Clearly, for all $x\in X,$ $T\left( x\right) \in X.$ Furthermore, 
\begin{equation*}
Karlstad=T\left( Karlstad\right)
\end{equation*}%
which is the unique fix-point of $T.$

Further more, we note that the distance between any two points $x,y$ in the
domain (large map), is strictly larger than the distance between $T\left(
x\right) $ and $T\left( y\right) ,$ i.e., between the corresponding points
in the small map.

Now, the contraction mapping theorem says that if;

\begin{enumerate}
\item For all $x\in X,T\left( x\right) \in X,$i.e., it maps elements in $X$
back into $X$\thinspace\ (Which in this case means that the small scale map
is completely inside the larger map), and

\item there is a distance measure $d,$such that for all $x,y$ $\in X,d\left(
x,y\right) >d\left( T\left( x\right) ,T\left( y\right) \right) $ (the small
scale map has a strictly smaller scale),
\end{enumerate}

then there is exactly one fix-point $x=T\left( x\right) .$

$T$ is then called a \emph{contraction mapping.}

Below, we generalize this result to another space $X$, namely the space of
continuous and bounded functions. Clearly, the Bellman $T$-operator maps
elements of the set bounded continuos (value)functions, which we call $C,$
into the same set.

What remains is to define a distance measure. We use the so called sup-norm.
This is defined as follows. Consider two functions (in our case value
functions), $w\left( s\right) $ and $v\left( s\right) ,s\in S.$ This two
functions maps elements in a state space $S\mathbf{\subseteq R}%
^{n}\rightarrow \mathbf{R.}$ $w$ and $v$ are themselves elements of $C.$
Then, the distance $d\left( w,v\right) $ between the two functions is the
maximum absolute difference, 
\begin{equation*}
d\left( w,v\right) \equiv \sup_{s\in \mathbf{S}}\left\Vert w\left( s\right)
-v\left( s\right) \right\Vert =\sup_{s\in \mathbf{S}}\sqrt{\left( w\left(
s\right) -v\left( s\right) \right) ^{2}}
\end{equation*}

Typically, our Bellman \thinspace $T$-operator is a contraction mapping if
we make sure we have bounded state spaces and there is strict discounting.
The next sections make these claims a little more formally.

\subsubsection{Complete Metric Spaces and Cauchy Sequences}

Let $\mathbf{X}$ be a vector space, i.e., a set on which addition and scalar
multiplication is defined. Also define an operator $d$: which we can think
of as measuring the (generalized) distance between any two elements of $%
\mathbf{X}$. We call $d$ a norm assumed to satisfy

\begin{enumerate}
\item Positivity $\forall x,y\,\in \mathbf{X}$, $d\left( x,y\right) \geq 0$
and $d\left( x,y\right) =0\Rightarrow x=y.$

\item Symmetry $\forall x,y\,\in \mathbf{X}$, $d\left( x,y\right) =d\left(
y,x\right) .$

\item Triangle inequality $\forall x,y,z\,\in \mathbf{X}$, $d\left(
x,z\right) \leq d\left( x,y\right) +d\left( y,z\right) $
\end{enumerate}

Now, we call $(\mathbf{X},d)$ a \emph{normed vector space} or a \emph{metric
space}. An example of such a space would be $\mathbf{R}^{n}$ together with
the Euclidean norm $d\left( x,y\right) \equiv \left\Vert x,y\right\Vert $.
Another example is the space $\mathbf{C(S)}$ of continuous, bounded
functions where each element is a function from $\mathbf{S\subseteq R}%
^{n}\rightarrow \mathbf{R}$ together with the \textquotedblright
sup-norm\textquotedblright\ defined as follows. For any two elements in $%
\mathbf{C}\left( \mathbf{S}\right) $, i.e., any two functions $w\left(
s\right) $ and $v\left( s\right) $, the distance $d$ between them is the
maximal euclidean distance, i.e., 
\begin{equation}
d\left( w,v\right) \equiv \sup_{s\in \mathbf{S}}\left\Vert w\left( s\right)
,v\left( s\right) \right\Vert
\end{equation}

Now let us define a Cauchy sequence. Intuitively, this is a sequence of
elements $\{x_{n}\}$ in a space $\mathbf{X}$ that come closer and closer to
each other, using some particular norm. More precisely, $\{x_{n}\}$ is
defined as a sequence of elements in $\mathbf{X}$ such that for all $%
\varepsilon >0$, there exist a number $n$, such that for all $m,p\geq n$, $%
d(x_{m},x_{p})<\varepsilon $. An example of such a sequence would be the
sequence $\{1,1/2,1/3,\}$ which is a Cauchy sequence using the Euclidean
norm. A Cauchy sequence converges if there is an element $y\in \mathbf{X}$
such that $\lim_{n\rightarrow \infty }d(x_{n},y)=0$. It may, of course, be
the case that the Cauchy sequence does not converge to a point in $\mathbf{X}
$. An example would be if we let $\mathbf{X}$ be the open interval $(0,1]$
and look at the Cauchy sequence $\{1,1/2,1/3,\}$ which is in $\mathbf{X}$
but converges to zero which is not in $\mathbf{X.}$

\subsubsection{Complete metric spaces}

Now we are ready to define the \emph{complete metric space}. This is a
metric space in which all Cauchy sequences converge to a point in the space.

\subsubsection{ Contraction Mapping}

Consider a metric space $(\mathbf{X},d)$ and look at an operator $T$ that
maps $\mathbf{X}$ $\rightarrow $ $\mathbf{X}$. $T$ is a contraction mapping
by definition if there exists a non-negative number $\rho \in \lbrack 0,1),$
such that for all elements $x,y\in \mathbf{X,}$

\begin{equation}
d\left( T\left( x\right) ,T\left( y\right) \right) \leq \rho d\left(
x,y\right) ,  \label{contractiondef}
\end{equation}
where we note that $\rho $ must be \emph{strictly} smaller than one.

\subsubsection{The Contraction Mapping Theorem}

Now we can state the very important \emph{contraction mapping theorem.}

\begin{result}
Consider a complete metric space $\left( \mathbf{X},d\right) $ and let $T:%
\mathbf{X\rightarrow X}$ be a contraction mapping. Then, $T$ has one unique
fixed point $x^{\ast }\in \mathbf{X}$, i.e., the solution to $x=T(x)$ always
exists and is unique. Furthermore, the sequence $x_{0}$, $T(x_{0})$,$%
T^{2}(x_{0})$,...,$T^{n}(x_{0})$ converges to $x^{\ast }$ for all $x_{0}\in 
\mathbf{X}$.
\end{result}

There are theorems that can be used to show that $T$ is a contraction
mapping.

\begin{result}
\label{Blackwell}Let the state space $\mathbf{S}$ be a subset of $\mathbf{R}%
^{n}$ and $\mathbf{B}(\mathbf{S})$ the set of all continuous, \emph{bounded}
functions from $\mathbf{S}$ to $\mathbf{R}$. Endowed with the sup-norm, $%
\mathbf{B}(\mathbf{S})$, is a complete metric space. Let $T$ be a map that
maps all elements of $\mathbf{B}(\mathbf{S)\rightarrow }$ $\mathbf{B}(%
\mathbf{S)}$. Then, $T$ is a contraction mapping if

\begin{enumerate}
\item for any functions $w(s),v(s)\in \mathbf{B}(\mathbf{S})$, the following
holds; if $w(s)\geq v(s)\forall s\in \mathbf{S}$ then $T\left( w(s)\right)
\geq T\left( v(s)\right) \forall s\in \mathbf{S}$ (\emph{monotonicity}), and

\item there is a $\beta \in \lbrack 0,1)$ such that for any constant $\kappa
\in R$, and any function $w\left( s\right) \in \mathbf{B}\left( \mathbf{S}%
\right) ,T\left( w(s)+\kappa \right) =T\left( w(s)\right) +\beta \kappa $. (%
\emph{discounting}).
\end{enumerate}
\end{result}

Usually it is straightforward to apply the previous result to show that if
we have strict discounting, the Bellman equation is a contraction mapping.
There is one major limitation which we have to live with, however, result %
\ref{Blackwell} and variants of it require bounded value functions.

Let us look at an example, where we apply result \ref{Blackwell}. Consider a
simple growth model, 
\begin{eqnarray}
&&\max_{\{c_{t}\}_{0}^{\infty }}E_{0}\sum_{t=0}^{\infty }\beta ^{t}u\left(
c_{t}\right) \\
\text{s.t. }k_{t+1} &=&f\left( k_{t}\right) -c_{t},\forall t\geq 0,k_{0}\,%
\text{given, and} \\
c_{t},k_{t} &>&0\text{ }\forall t\text{.}
\end{eqnarray}
where $u$ is a continuous and increasing (utility) function with $u\left(
0\right) \geq u_{\min }>-\infty $ and $0\leq \beta <1.$ To use the theorems
we need to make some assumptions. First, we need boundedness. For this
purpose, we assume 
\begin{eqnarray}
f\left( 0\right) &=&0, \\
f^{\prime }\left( k\right) &\geq &0\forall k\geq 0, \\
\exists \bar{k} &>&0,\text{such that }f\left( k\right) \leq k\forall k\geq 
\bar{k}.
\end{eqnarray}

Now, define $\mathbf{S}=[0,\bar{k}]$ and note that if $k_{0}$ is in $\mathbf{%
S}$, so is all admissible $k_{t}$. Then, $u\left( c_{t}\right) $ is bounded
since $u_{\min }\leq u\left( c\right) \leq u\left( \bar{k}\right) $ implying
that any value function must satisfy $\frac{u_{\min }}{1-\beta }\leq
V(k)\leq \frac{u\left( \bar{k}\right) }{1-\beta }.$ By restricting the state
space $\mathbf{S}=[0,\bar{k}]$, we can therefore restrict our search for
value functions that are bounded on our state space.

Now, let us establish that the following mapping is a contraction: 
\begin{equation}
V\left( k\right) =\max_{c\geq 0}u\left( c\right) +\beta V\left( f\left(
k)-c\right) \right) \equiv T\left( V\left( k\right) \right)
\end{equation}

Regarding condition 1, we need for any two bounded functions $v\left(
k\right) ,w\left( k\right) ,k\in S,.$ 
\begin{equation}
v\left( k\right) \geq w\left( k\right) \forall k\Rightarrow T\left( v\left(
k\right) \right) \geq T\left( w\left( k\right) \right) \forall k.
\end{equation}%
which is satisfied. To see this, define 
\begin{equation}
c^{\ast }=\arg \max_{c}u\left( c\right) +\beta w\left( f\left( k\right)
-c\right)
\end{equation}%
then, 
\begin{eqnarray}
T\left( v\left( k\right) \right) &=&\max_{c}u\left( c\right) +\beta v\left(
f\left( k\right) -c\right) \geq u\left( c^{\ast }\right) +\beta v\left(
f\left( k\right) -c^{\ast }\right) \\
&\geq &u\left( c^{\ast }\right) +\beta w\left( f\left( k\right) -c^{\ast
}\right) =T\left( w\left( k\right) \right) .
\end{eqnarray}

Regarding condition 2, we note 
\begin{eqnarray}
T\left( v\left( k\right) +\kappa \right) &=&\max_{c\geq 0}u\left( c\right)
+\beta \left( v\left( f\left( k-c\right) \right) +\kappa \right) \\
&=&\max_{c\geq 0}u\left( c\right) +\beta v\left( f\left( k-c\right) \right)
+\beta \kappa \\
&=&T\left( v\left( k\right) \right) +\beta \kappa ,
\end{eqnarray}%
so the second condition is satisfied too. Thus, the Bellman equation is a
contraction mapping and always has one and one only solution, $V\left(
k\right) .$

\newpage

\section{Dynamic Optimization in Continuous Time}

\subsection{Dynamic programming in continuous time}

%TCIMACRO{\TeXButton{TeX field}{\setcounter{equation}{0}}}%
%BeginExpansion
\setcounter{equation}{0}%
%EndExpansion
Consider the problem 
\begin{eqnarray}
&&\max_{k\left( t\right) _{0}^{T},c\left( t\right)
_{0}^{T}}\int_{0}^{T}e^{-rt}u\left( k,c,t\right) dt  \label{Conttime1} \\
s.t.\,\dot{k} &=&f\left( k,c,t\right)  \label{LoM1} \\
k\left( 0\right) &=&\underline{k},
\end{eqnarray}

with 
\begin{gather}
k\left( T\right) =\bar{k}\text{ (case 1), or}  \label{endcases} \\
\text{ }k\left( T\right) \text{ free (case 2), or} \\
k\left( T\right) \geq \bar{k}\text{ (case 3).}
\end{gather}

Thinking of the integral in the maximand as a sum of rectangles with base $%
dt $ and height $e^{-rt}u\left( k,c,t\right) $, we can approximate the
problem with a discrete time problem. Noting that for a small time interval $%
dt$, $k\left( t+dt\right) =k\left( t\right) +f\left( k,c,t\right) dt,$ we
can write the current value Bellman equation 
\begin{equation}
V\left( k,t\right) =\max_{c}\left\{ u\left( k,c,t\right) dt+e^{-rdt}V\left(
k+f\left( k,c,t\right) dt,t+dt\right) \right\} .
\end{equation}

We can then make Taylor approximations; 
\begin{gather}
V\left( k,t\right) =\max_{c}{\Large \{}u\left( k,c,t\right) dt \\
+\left( 1-rdt\right) \left( V\left( k,t\right) +V_{k}\left( k,t\right)
f\left( k,c,t\right) dt+V_{t}\left( k,t\right) dt\right) {\Large \}}
\end{gather}

Subtracting $\left( 1-rdt\right) V\left( k,t\right) $ from both sides,
dividing by $dt$ and then letting $dt$ go to zero yields, 
\begin{eqnarray}
rV\left( k,t\right) &=&\max_{c}\left\{ u\left( k,c,t\right) +V_{k}\left(
k,t\right) f\left( k,c,t\right) +V_{t}\left( k,t\right) \right\}
\label{arbitrage} \\
&=&\max_{c}\left\{ u\left( k,c,t\right) +\frac{dV\left( k,t\right) }{dt}%
\right\} ,  \label{arbitrageb}
\end{eqnarray}%
where we note that $dV\left( k,t\right) /dt$ is the total time derivative of
the value function. Sometimes, this equation is referred to as an asset
pricing equation -- and interpreted as follows; if an asset is correctly
valued (not providing arbitrage opportunities) the opportunity cost of
holding it (the LHS of (\ref{arbitrage})) equals the (optimal) sum of the
immediate pay-off or dividend (the first term of (\ref{arbitrageb})) and the
capital gain (the second term of (\ref{arbitrageb})).

Sometimes we can use the guess and verify technique to solve for the value
function if the problem is autonomous. An alternative is to use \emph{%
Optimal Control} and \emph{Pontryagin's maximum principle.}

\subsection{Optimal Control}

Consider the problem in (\ref{Conttime1}) with the continuous time Bellman
equation (\ref{arbitrage}). Noting that since the term $V_{t}\left(
k,t\right) $ is independent of $c,$ we can rewrite (\ref{arbitrage}) 
\begin{equation}
rV\left( k,t\right) -V_{t}\left( k,t\right) =\max_{c}\left\{ u\left(
k,c,t\right) +V_{k}\left( k,t\right) f\left( k,c,t\right) \right\}
\label{Arbitrage2}
\end{equation}

Defining the \emph{co-state} or current shadow value variable as the
derivative of the value function w.r.t. $k,$ along its optimal path $k^{\ast
}$, 
\begin{equation}
\lambda \left( t\right) \equiv V_{k}\left( k^{\ast },t\right) .
\end{equation}

We see that a necessary condition for $c^{\ast }\left( t\right) $ to be
optimal is that it is given by 
\begin{equation}
c^{\ast }\left( t\right) =\arg \max_{c}\left\{ u\left( k,c,t\right) +\lambda
\left( t\right) f\left( k,c,t\right) \right\} .  \label{Pont1b}
\end{equation}

Now, we need to pin down $\lambda \left( t\right) .$ For this purpose we
analyze how $\lambda $ develops over time by deriving a differential
equation for $\lambda \left( t\right) .$Taking derivatives w.r.t. $k$ of the
identity (\ref{Arbitrage2}), we get 
\begin{gather}
rV_{k}\left( k,t\right) -V_{kt}\left( k,t\right) =u_{k}\left( k,c^{\ast
},t\right) +V_{kk}\left( k,t\right) f\left( k,c^{\ast },t\right) \\
+V_{k}\left( k,t\right) f_{k}\left( k,c^{\ast },t\right) , \\
r\lambda \left( t\right) -\left( V_{tk}\left( k,t\right) +V_{kk}\left(
k,t\right) \dot{k}\right) =u_{k}\left( k,c^{\ast },t\right) +\lambda \left(
t\right) f_{k}\left( k,c^{\ast },t\right) , \\
r\lambda \left( t\right) -\dot{\lambda}\left( t\right) =u_{k}\left(
k,c^{\ast },t\right) +\lambda \left( t\right) f_{k}\left( k,c^{\ast
},t\right) ,  \label{Pont2}
\end{gather}%
where in second equation, we changed the order of diffentiation, using the
fact that $\frac{\partial }{\partial k}\frac{\partial V\left( k,t\right) }{%
\partial t}=\frac{\partial }{\partial t}\frac{\partial V\left( k,t\right) }{%
\partial k}.$

Equation (\ref{Pont1b}) and (\ref{Pont2}) form the basis of Pontryagin's
maximum principle. According to this, we derive necessary (and sometimes
sufficient) conditions for an optimal control. We do this without explicitly
solving for the value function. We first define the \emph{current value
Hamiltoninan}. This is the sum of the instantaneous payoff and the product
of costate(s) and the function determining the law-of-motion of the state
variable; 
\begin{equation}
\mathcal{H}\left( k,c,\lambda ,t\right) \mathcal{\equiv }u\left(
k,c,t\right) +\lambda \left( t\right) f\left( k,c,t\right) .
\end{equation}

Note that the Hamiltonian has the same interpretation as the RHS\ of the
Bellman equation without the max-operator. In words, it is the sum of the
flow of current pay-off and the generation of future pay-off. In the Bellman
equation, we use the value function to measure future payoffs while the
Hamiltonian uses the shadow-value $\lambda .$

According to Pontryagin's maximum principle, the optimal control $c^{\ast
}\left( t\right) $ maximizes the Hamiltonian at each instant, the co-state
(or shadow value) satisfies the differential equation, $r\lambda \left(
t\right) -\dot{\lambda}\left( t\right) =\mathcal{H}_{k}\left( k,c,\lambda
,t\right) $. These necessary conditions will provide differential equations
which we need to solve. Typically we have one initial condition for each
state variable. But we need more information to solve the system since we
also have the control variable(s). In case (1) of (\ref{endcases}), we have
the necessary additional information. In case (2), it must be that the
shadow value of the state variable approach zero as $t\rightarrow T.$ This
is called, \emph{the transversality condition. }In case (3), either the
inequality is slack, in which case $\lambda \left( T\right) =0,$ or it
binds, giving the necessary additional info in both cases. We can the
summarize:

\begin{result}
According to the \emph{Pontryagin's maximum principle}\newline
1. An optimal control $c^{\ast }\left( t\right) ,$ satisfies 
\begin{equation}
c^{\ast }\left( t\right) =\arg \max_{c}\mathcal{H}\left( k,c,\lambda
,t\right) ,  \label{pontr1}
\end{equation}%
\newline
2. where the current co-state or shadow value $\lambda \left( t\right) $ ,
is a continuous function of time, satisfying the differential equation. 
\begin{equation}
r\lambda \left( t\right) -\dot{\lambda}\left( t\right) =\mathcal{H}%
_{k}\left( k,c,\lambda ,t\right) ,  \label{pont2}
\end{equation}%
except at points in time where $c$ is discontinuous and with \newline
3. end-condition(s) provided by 
\begin{align}
k\left( T\right) & =\bar{k}\text{ (case 1) or}  \label{pont3a} \\
\lambda \left( T\right) & =0\text{ if }k\left( T\right) \text{ is free (case
2), or}  \label{pont3b} \\
\lambda \left( T\right) & \geq 0,\text{ and }\lambda \left( T\right) \left(
k\left( T\right) -\bar{k}\right) =0\text{ (case (3).}  \label{pont3c}
\end{align}%
\newline
In addition, the path for $k\left( t\right) $ must satisfy the initial
condition, $k\left( 0\right) =\underline{k}$, and $\dot{k}\left( t\right) =%
\mathcal{H}_{\lambda }\left( k,c,\lambda ,t\right) =f\left( k,c,t\right) .$
\end{result}

\subsubsection{The consumption problem}

As an example, consider the problem standard consumption-savings (Ramsey)
problem 
\begin{eqnarray}
&&\max_{k\left( t\right) _{0}^{T},c\left( t\right)
_{0}^{T}}\int_{0}^{T}e^{-rt}u\left( c\right) dt  \label{Ramseyob} \\
s.t.\,\dot{k} &=&f\left( k\right) -c \\
k\left( 0\right) &=&\underline{k}, \\
k\left( T\right) &\geq &0,
\end{eqnarray}%
with $u$ and $f$ increasing and concave. The current value Hamiltonian is 
\begin{equation}
\mathcal{H}\left( k,c,\lambda ,t\right) =u\left( c\right) +\lambda \left(
t\right) \left( f\left( k\right) -c\right) .
\end{equation}

Thus, $c^{\ast }\left( t\right) =\arg \max_{c}u\left( c\right) +\lambda
\left( t\right) \left( f\left( k\right) -c\right) $ which is interior,
implying 
\begin{gather}
\mathcal{H}_{c}\left( k,c,\lambda ,t\right) =0. \\
\rightarrow u^{\prime }\left( c^{\ast }\right) =\lambda .
\end{gather}

Furthermore, from the second condition 
\begin{equation}
\mathcal{H}_{k}\left( k,c,\lambda ,t\right) =\lambda f^{\prime }\left(
k\right) =r\lambda -\dot{\lambda}.  \label{lam}
\end{equation}

Taking time-derivatives of $u^{\prime }\left( c^{\ast }\right) =\lambda $
and substituting into (\ref{lam}) we get 
\begin{eqnarray}
u^{\prime \prime }\left( c^{\ast }\right) \dot{c}^{\ast } &=&\dot{\lambda},
\\
u^{\prime }\left( c^{\ast }\right) f^{\prime }\left( k\right) &=&ru^{\prime
}\left( c^{\ast }\right) -u^{\prime \prime }\left( c^{\ast }\right) \dot{c}%
^{\ast }, \\
&\rightarrow &\dot{c}^{\ast }=\frac{u^{\prime }\left( c^{\ast }\right) }{%
-u^{\prime \prime }\left( c^{\ast }\right) }\left( f^{\prime }\left(
k\right) -r\right) ,
\end{eqnarray}%
which is the Euler equation we have seen before. To analyze the behavior of
this system, we can use the phase-diagram as in section \ref{phasediagram}.

To get a closed form solution, i.e., an expression for the endogenous
variables in terms of only the exogenous ones, we must specify the utility
and production functions. Considering first the utility functions, we have
two important special cases. First, CARA utility, 
\begin{eqnarray}
u\left( c\right) &=&-\frac{e^{-\gamma c}}{\gamma }, \\
&\rightarrow &u^{\prime }\left( c\right) =e^{-\gamma c},u^{\prime \prime
}\left( c\right) =-\gamma e^{-\gamma c}
\end{eqnarray}%
in which case we get 
\begin{equation}
\dot{c}^{\ast }=\frac{1}{\gamma }\left( f^{\prime }\left( k\right) -r\right)
,
\end{equation}%
i.e., consumption increase is a linear function in the difference between
the marginal return on savings and the subjective discount rate. The other
case is CRRA, 
\begin{eqnarray}
u\left( c\right) &=&\frac{c^{1-\sigma }}{^{1-\sigma }}, \\
&\rightarrow &u^{\prime }\left( c\right) =c^{-\sigma },u^{\prime \prime
}\left( c\right) =-\sigma c^{-\sigma -1},
\end{eqnarray}%
yielding 
\begin{eqnarray}
\dot{c}^{\ast } &=&\frac{c^{-\sigma }}{\sigma c^{-\sigma -1}}\left(
f^{\prime }\left( k\right) -r\right) =\frac{c}{\sigma }\left( f^{\prime
}\left( k\right) -r\right)  \label{Ramseyc} \\
\frac{\dot{c}^{\ast }}{c} &=&\frac{1}{\sigma }\left( f^{\prime }\left(
k\right) -r\right) ,  \label{Ramseyk}
\end{eqnarray}%
i.e., the consumption growth \emph{rate }is a linear function in $f^{\prime
}\left( k\right) -r.$ The sensitivity is given by $1/\sigma $, which we call
the intertemporal elasticity of substitution. Here, we also have that $%
\sigma $ is the constant of relative risk aversion.

What about the transversality condition? In this case, we know $u^{\prime
}\left( c^{\ast }\left( T\right) \right) =\lambda \left( T\right) .$ So
since utility in the examples is unbounded, i.e., $u^{\prime }\left(
c\right) >0$ for all finite $c,$ $\lambda \left( T\right) $ cannot be 0,
instead $k\left( T\right) $ is zero. In other words, whenever consumption is
valuable at $T$, the lower bound on $k$ should bind and nothing should be
left.

Let us complete the example by assuming, for simplicity, a linear (Romer
type) production function $f\left( k\right) =Ak$. In the CRRA\ case, we get
the linear system 
\begin{eqnarray}
\dot{c}^{\ast } &=&\frac{A-r}{\sigma }c \\
\dot{k} &=&-c+Ak \\
\left[ 
\begin{array}{c}
\dot{c}^{\ast } \\ 
\dot{k}%
\end{array}%
\right] &=&\left[ 
\begin{array}{cc}
\frac{A-r}{\sigma } & 0 \\ 
-1 & A%
\end{array}%
\right] \left[ 
\begin{array}{c}
c^{\ast } \\ 
k%
\end{array}%
\right] .
\end{eqnarray}

This system has roots $A$ and $\frac{A-r}{\sigma }$ and the matrix of
eigenvectors is 
\begin{equation}
\left[ 
\begin{array}{cc}
0 & \frac{r+A\left( \sigma -1\right) }{\sigma } \\ 
1 & 1%
\end{array}
\right] \equiv \mathbf{B}^{-1}.
\end{equation}

Consequently, the solution is of the form 
\begin{eqnarray}
\left[ 
\begin{array}{c}
c^{\ast }\left( t\right) \\ 
k\left( t\right)%
\end{array}%
\right] &=&\mathbf{B}^{-1}\left[ 
\begin{array}{cc}
e^{At} & 0 \\ 
0 & e^{\frac{A-r}{\sigma }t}%
\end{array}%
\right] \left[ 
\begin{array}{c}
\kappa _{1} \\ 
\kappa _{2}%
\end{array}%
\right] \\
&=&\left[ 
\begin{array}{cc}
0 & \frac{r+A\left( \sigma -1\right) }{\sigma }e^{\frac{A-r}{\sigma }t} \\ 
e^{At} & e^{\frac{A-r}{\sigma }t}%
\end{array}%
\right] \left[ 
\begin{array}{c}
\kappa _{1} \\ 
\kappa _{2}%
\end{array}%
\right] ,
\end{eqnarray}%
where $\kappa _{1}$ and $\kappa _{2}$ are two integration constants. We
solve for the latter by using $k\left( 0\right) =\underline{k}$ and $k\left(
T\right) =0.$%
\begin{eqnarray}
k\left( 0\right) &=&\underline{k}=\left[ 
\begin{array}{cc}
1 & 1%
\end{array}%
\right] \left[ 
\begin{array}{c}
\kappa _{1} \\ 
\kappa _{2}%
\end{array}%
\right] \\
&=&\kappa _{1}+\kappa _{2} \\
k\left( T\right) &=&0=\left[ 
\begin{array}{cc}
e^{T} & e^{\frac{A-r}{\sigma }T}%
\end{array}%
\right] \left[ 
\begin{array}{c}
\kappa _{1} \\ 
\kappa _{2}%
\end{array}%
\right] \\
&\rightarrow &\kappa _{1}=\frac{e^{\frac{A-r-\sigma }{\sigma }T}}{e^{\frac{%
A-r-\sigma }{\sigma }T}-1}\underline{k},\kappa _{2}=\frac{1}{1-e^{\frac{%
A-r-\sigma }{\sigma }T}}\underline{k}.
\end{eqnarray}

We can now, for example, evaluate 
\begin{eqnarray}
c^{\ast }\left( t\right) &=&\frac{r+A\left( \sigma -1\right) }{\sigma }e^{%
\frac{A-r}{\sigma }t}\frac{1}{1-e^{\frac{A-r-\sigma }{\sigma }T}}\underline{k%
}, \\
c^{\ast }\left( 0\right) &=&\frac{r+A\left( \sigma -1\right) }{\sigma }\frac{%
1}{1-e^{\frac{A-r-\sigma }{\sigma }T}}\underline{k}.
\end{eqnarray}%
$\allowbreak $

\subsection{Sufficiency}

Assume that $f$ and $u$ are concave in $k,c$ and $\lambda \geq 0$. This
implies that the Hamiltonian is concave in $k,c$. Then, Pontryagin's
necessary conditions (\ref{pontr1}) and (\ref{pont2}) and (\ref{pont3a}), or
(\ref{pont3b}) or (\ref{pont3c}) are sufficient.

\subsection{Infinite horizon}

Consider the infinite horizon problem 
\begin{eqnarray}
&&\max_{c\left( t\right) _{0}^{\infty }}\int_{0}^{\infty }e^{-rt}u\left(
k,c,t\right) dt \\
s.t.\,\dot{k} &=&f\left( k,c,t\right) \\
k\left( 0\right) &=&\underline{k},
\end{eqnarray}

Pontryagin's conditions (\ref{pontr1}) and (\ref{pont2}) are necessary also
in the infinite horizon case, provided, of course, that there is a well
defined solution. If there is a binding restriction on the state variable of
the type $\lim_{T\rightarrow \infty }k\left( T\right) =\bar{k}$, this can
help us pin down the solution. The finite horizon transversality conditions
can, however, not immediately be used in the infinite horizon case. Suppose
the \emph{maximized} Hamiltonian is concave in $k$ for every $t$, then the
conditions (\ref{pontr1}) and (\ref{pont2}) plus the infinite horizon
transversality conditions 
\begin{eqnarray}
\lim_{T\rightarrow \infty }e^{-rT}\lambda \left( T\right) k\left( T\right)
&=&0,\text{ and }  \label{trans1} \\
\lim_{T\rightarrow \infty }e^{-rT}\lambda \left( T\right) &\geq &0,
\label{trans2}
\end{eqnarray}%
provide a \emph{sufficient} set of conditions for optimality (see de la
Fuente, 1999, p 577). Often, the Hamiltonian is concave in $k,c$ together.
This is sufficient for the maximized Hamiltonian to be concave in $k.$

Sometimes a so called \emph{No-Ponzi} condition helps us to make sure that
the transversality conditions are satisfied. Suppose, for example, the
pay-off $u\left( k,c,t\right) =u\left( c\right) $, that $k$ represents debt
of the agent and for simplicity that $f\left( k,c,t\right) =c+\rho k-w$, so
debt increases by the difference between consumption plus interest payments $%
\rho k$ and the wage $w$. It is reasonable to assume that creditors demand
to be repaid in a present value sense -- the discounted value of future
repayment should always be at least as large as debt. This is the No-Ponzi
condition. When in addition, the agent prefers to pay back no more than he
ows, the implication is. 
\begin{equation}
\lim_{T\rightarrow \infty }e^{-\rho T}k\left( T\right) =0.
\end{equation}%
To see this, solve 
\begin{equation}
\dot{k}\left( t\right) -\rho k\left( t\right) =c\left( t\right) -w\left(
t\right)
\end{equation}%
giving 
\begin{eqnarray}
e^{-\rho t}\left( \dot{k}\left( t\right) -\rho k\left( t\right) \right)
&=&e^{-\rho t}\left( c\left( t\right) -w\left( t\right) \right) \\
e^{-\rho t}k\left( t\right) &=&\int_{0}^{t}\left( e^{-\rho s}\left( c\left(
t\right) -w\left( t\right) \right) \right) ds+k\left( 0\right) \\
\lim_{T\rightarrow \infty }e^{-\rho T}k\left( T\right) &=&\int_{0}^{\infty
}\left( e^{-\rho s}\left( c\left( t\right) -w\left( t\right) \right) \right)
ds+k\left( 0\right)
\end{eqnarray}%
So, the No-Ponzi requirement is that if the PDV of "mortgage" repayments is
no smaller than initial debt, i.e., $-\int_{0}^{\infty }\left( e^{-\rho
s}\left( c\left( t\right) -w\left( t\right) \right) \right) ds$ $k\left(
0\right) \geq k\left( 0\right) ,$ then $\lim_{T\rightarrow \infty }e^{-\rho
T}k\left( T\right) \leq 0.$ Clearly, when marginal utility is strictly
positive, the individual would never want to satisfy this with inequality,
since he could then increase consumption. Therefore, $\lim_{T\rightarrow
\infty }e^{-\rho T}k\left( T\right) =0.$

The second necessary condition in (\ref{pont2}) is now 
\begin{eqnarray}
-\left( \rho -r\right) \lambda \left( t\right) &=&\dot{\lambda}\left(
t\right) \\
\lambda \left( t\right) &=&\lambda \left( 0\right) e^{-\left( \rho -r\right)
t}.
\end{eqnarray}%
So provided marginal utility is positive at $t=0$, (\ref{trans2}) is
satisfied. Furthermore, 
\begin{eqnarray}
\lim_{T\rightarrow \infty }e^{-rT}\lambda \left( T\right) k\left( T\right)
&=&\lambda \left( 0\right) \lim_{T\rightarrow \infty }e^{-rT}e^{-\left( \rho
-r\right) T}k\left( T\right) \\
&=&\lambda \left( 0\right) \lim_{T\rightarrow \infty }e^{-\rho T}k\left(
T\right) \\
&=&0.
\end{eqnarray}%
where the last equality is the No-Ponzi condition.

Sometimes, the sufficient conditions allow us to identify the optimal
control as the stable manifold (saddle-path) leading to a saddle-point
stable steady state. Consider again the problem 
\begin{eqnarray}
&&\max_{c\left( t\right) _{0}^{\infty }}\int_{0}^{\infty }e^{-rt}u\left(
c\right) dt \\
\text{s.t.}\,\dot{k} &=&f\left( k\right) -c,\text{ }k\left( 0\right) =%
\underline{k},
\end{eqnarray}%
which we graphically analyzed in section \ref{phasediagram} showing the
existence of saddle-path and a steady state with 
\begin{eqnarray}
f^{\prime }\left( k^{ss}\right) &=&r \\
c^{ss} &=&f\left( k^{ss}\right)
\end{eqnarray}

Restating the current value Hamiltonian 
\begin{equation}
\mathcal{H}\left( k,c,\lambda ,t\right) =u\left( c\right) +\lambda \left(
t\right) \left( f\left( k\right) -c\right) ,
\end{equation}
we note that if both $u\left( c\right) $ and $f\left( k\right) $ are
concave, and $\lambda \left( t\right) \geq 0,$ $\mathcal{H}\left(
k,c,\lambda ,t\right) $ is concave in $k,c$ so the conditions for using the
sufficiency result are satisfied. In addition to (\ref{Ramseyc}), (\ref%
{Ramseyk}), and $k\left( 0\right) =\underline{k}$, we thus only need to
verify that (\ref{trans1}) and (\ref{trans2}) are satisfied. This is
straightforward, 
\begin{eqnarray}
\lim_{T\rightarrow \infty }e^{-rT}\lambda \left( T\right)
&=&\lim_{T\rightarrow \infty }e^{-rT}u^{\prime }\left( c^{ss}\right) =0, \\
\lim_{T\rightarrow \infty }e^{-rT}\lambda \left( T\right) k\left( T\right)
&=&\lim_{T\rightarrow \infty }e^{-rT}u^{\prime }\left( c^{ss}\right)
k^{ss}=0.
\end{eqnarray}

\subsection{Present value Hamiltonian}

Sometimes, it is convenient to define the present value Hamiltonian, i.e.,
expressing everything in values as seen from time 0. In problem (\ref%
{Conttime1}), the present value Hamiltonian is given by 
\begin{equation}
\mathbf{H}\left( k,c,\mu ,t\right) =e^{-rt}u\left( k,c,t\right) +\mu f\left(
k,c,t\right) ,
\end{equation}
where $\mu \left( t\right) $ is the present shadow value of the state
variable. In this case, the necessary conditions for optimality are 
\begin{eqnarray}
c^{\ast }\left( t\right) &=&\arg \max_{c}\mathbf{H}\left( k,c,\mu ,t\right)
\\
-\dot{\mu}\left( t\right) &=&\mathbf{H}_{k}\left( k,c,\mu ,t\right) \\
\dot{k} &=&\mathbf{H}_{\mu }\left( k,c,\mu ,t\right)
\end{eqnarray}

In the finite horizon case, the transversality conditions are the same in
terms of $\mu \left( T\right) $ and $\lambda \left( T\right) .$In the
infinite horizon case, we note that 
\begin{equation}
\mu \left( T\right) =e^{-rT}\lambda \left( T\right) ,
\end{equation}%
so the conditions (\ref{trans1}) and (\ref{trans2}) become 
\begin{eqnarray}
\lim_{T\rightarrow \infty }\mu \left( T\right) k\left( T\right) &=&0,\text{
and } \\
\lim_{T\rightarrow \infty }\mu \left( T\right) &\geq &0.
\end{eqnarray}

Note also that Michel (Econometrica, July 1982) states 
\begin{equation*}
\lim_{\lim_{T\rightarrow \infty }}\mathbf{H}\left( k,c,\mu ,t\right) =0
\end{equation*}%
as a necessary general transversality condition in the infinite horizon.
This condition seems intuitive in the sense that it can be interpreted as
saying that all profit opportunities should be exhausted in the long run.
However, the issue of general transversaility conditions is still in
dispute, and caution is warranted although counterexamples typically seem
quite "cooked".

\subsection{Bang-Bang solutions}

We have already noted that the condition $r\lambda \left( t\right) -\dot{%
\lambda}\left( t\right) =\mathcal{H}_{k}\left( k,c,\lambda ,t\right) $ and $-%
\dot{\mu}\left( t\right) =\mathbf{H}_{k}\left( k,c,\mu ,t\right) $ holds
everywhere except at points in time where $c$ is discontinuous. Let us take
an example when such discontinuities arise because the Hamiltonian is linear
in the control. The control variable has to be chosen from a compact set and
the linearity implies that the optimal control jumps from one end of the
constraint to the other -- a \emph{bang-bang} solution.

An individual allocates his daily time between working and enjoying leisure.
When enjoying leisure, his utility increases in the amount of toys he has to
play with. If he works, he gets money which he spends on increasing his
stock of toys. On the other hand, while working, he cannot enjoy his toys.
The problem can be written, 
\begin{eqnarray*}
&&\max_{\{l\}_{0}^{T}}\int\limits_{0}^{T}\left( 1-l\left( t\right) \right)
U\left( x\left( t\right) \right) dt \\
s.t.\text{ }\dot{x}(t) &=&l(t), \\
1-l(t) &\geq &0, \\
l(t) &\geq &0, \\
x(0) &=&0,
\end{eqnarray*}%
where $l\in \lbrack 0,1]$ is the share of time per day spent on working, $x$
is the stock of toys and $U$ is a strictly increasing utility function with $%
U(0)=0$.

The Hamiltonian for the problem is%
\begin{equation*}
H=\left( 1-l\left( t\right) \right) U\left( x\left( t\right) \right)
+\lambda \left( t\right) \left( l(t)\right) .
\end{equation*}

Since there is no terminal condition on $x\left( T\right) ,$ the
transversality condition is $\lambda \left( T\right) =0.$The necessary
conditions for optimality are%
\begin{eqnarray}
l &=&\arg \max_{l}\left( 1-l\left( t\right) \right) U\left( x\left( t\right)
\right) +\lambda \left( t\right) \left( l(t)\right) , \\
H_{x} &=&\left( 1-l\left( t\right) \right) U_{x}\left( x\left( t\right)
\right) =-\dot{\lambda}\left( t\right) .  \label{eq_Bang2}
\end{eqnarray}

Clearly, $\lambda >U\left( x\right) $ implies $l^{\ast }=1$ and vice versa.
From (\ref{eq_Bang2}) we see that $\lambda $ is weakly falling. Therefore, $%
\lambda \left( 0\right) >0$ (Otherwise, the total payoff would be zero which
cannot be optimal). Since $U\left( x\left( 0\right) \right) =0,$ $l\left(
0\right) =1.$

Furthermore, it cannot be optimal to work all time. We conclude that the
optimal policy is a stopping rule -- work, i.e., set $l\left( t\right) =1$
until $t^{\ast }$ then set $l\left( t\right) =0$ for $t\in \left[ t^{\ast },T%
\right] .$ From (\ref{eq_Bang2}) we see that $\dot{\lambda}\left( t\right)
=0 $ for $t\in \lbrack 0,t^{\ast })$ and we also have $\dot{x}(t)=1$ for $%
t\in \lbrack 0,t^{\ast }).$ For $t\in (t^{\ast },T],-\dot{\lambda}\left(
t\right) =-\frac{1}{2}x\left( t^{\ast }\right) ^{-1/2}$ with end condition $%
\lambda \left( T\right) =0.$ This in combination with the fact that $\lambda 
$ is continuos and $\lambda \left( t^{\ast }\right) =U\left( t^{\ast
}\right) $ is enough to solve the problem.

If we specify for example $U(x)=\sqrt{x},$ we get a closed form solution.
Summarizing 
\begin{eqnarray*}
x\left( t\right) &=&\left\{ 
\begin{array}{c}
t\text{ for }t<t^{\ast } \\ 
t^{\ast }\text{ for }t>t^{\ast }%
\end{array}%
\right. \\
x\left( 0\right) &=&0 \\
\dot{\lambda}\left( t\right) &=&\left\{ 
\begin{array}{c}
0\text{ for }t<t^{\ast } \\ 
-U_{x}\left( x\left( t\right) \right) =-\frac{1}{2}x\left( t^{\ast }\right)
^{-1/2}\text{for }t>t^{\ast }%
\end{array}%
\right. \\
\lambda \left( t\right) &=&\left\{ 
\begin{array}{c}
\lambda \left( t^{\ast }\right) \text{ for }t<t^{\ast } \\ 
\text{ }\lambda \left( t^{\ast }\right) -\frac{1}{2}x\left( t^{\ast }\right)
^{-1/2}\left( t-t^{\ast }\right) \text{ for }t>t^{\ast }%
\end{array}%
\right. \\
\lambda \left( t^{\ast }\right) &=&U\left( x\left( t^{\ast }\right) \right)
=\left( t^{\ast }\right) ^{\frac{1}{2}} \\
\lambda \left( T\right) &=&0.
\end{eqnarray*}

Using $\lambda \left( t^{\ast }\right) =U\left( x\left( t^{\ast }\right)
\right) =\left( t^{\ast }\right) ^{\frac{1}{2}}$ in $\lambda \left( t^{\ast
}\right) -\frac{1}{2}x\left( t^{\ast }\right) ^{-1/2}\left( t-t^{\ast
}\right) $ and solving for $t^{\ast }$ gives, 
\begin{equation*}
t^{\ast }=\frac{1}{3}T.
\end{equation*}

\subsection{Many state variables and controls}

Having several state variables and controls pose no principle problem.
Neither does pointwise discontinuities in the control variable. To
generalize, suppose we have $n$ state variables and $n$ controls. An optimal
control maximizes the Hamiltonian over all available controls $\mathbf{c}$%
\begin{equation}
\mathbf{c}^{\ast }\left( t\right) =\arg \max_{\mathbf{c}}\mathcal{H}\left( 
\mathbf{k,c},\mathbf{\lambda },t\right) \equiv u\left( \mathbf{k,c},t\right)
+\sum_{i=1}^{n}\lambda _{i}f_{i}\left( \mathbf{k,c},t\right) .
\end{equation}%
where $\lambda _{i}$ is the shadow value associated with the state variable $%
k_{i}$. Each $\lambda _{i}\left( t\right) $ is continuous and satisfies the
differential equation 
\begin{equation}
r\lambda _{i}\left( t\right) -\dot{\lambda}_{i}\left( t\right) =\frac{%
\partial }{\partial k_{i}}\mathcal{H}\left( \mathbf{k,c},\mathbf{\lambda }%
,t\right) ,
\end{equation}%
except when $c$ is discontinuous. For the transversality conditions, we have 
\begin{align}
k_{i}\left( T\right) & =\bar{k}_{i}\text{ case (1), or} \\
\lambda _{i}\left( T\right) & =0\text{ case (2), or} \\
\lambda _{i}\left( T\right) & \geq 0,\text{ and }\lambda _{i}\left( T\right)
\left( k_{i}\left( T\right) -\bar{k}_{i}\right) =0\text{ case (3),}
\end{align}%
for end-conditions for state variable $i$ belonging to case 1, 2 or 3.

\newpage

\section{Some numerical methods}

%TCIMACRO{\TeXButton{TeX field}{\setcounter{equation}{0}}}%
%BeginExpansion
\setcounter{equation}{0}%
%EndExpansion

\subsection{Numerical solution to Bellman equations}

When we cannot solve the Bellman equation analytically, there are several
methods to approximate a solution numerically. One of the most
straightforward methods when the problem is autonomous, is to discretize the
state space and iterate on the Bellman equation until it converges. When the
Bellman equation is a contraction mapping, strong results make sure that
this procedure converges to the correct value function.

When we discretize the state space, we restrict the state variable to take
values from a finite set 
\begin{equation}
k_{t}\in \left[ k^{1},k^{2},...k^{n}\right] \equiv K
\end{equation}%
where the superscripts index the elements of $K.$ We then solve for $c$ from
the law-of-motion for $k$%
\begin{eqnarray}
k_{t+1} &=&f\left( k_{t},c_{t}\right) \\
&\implies &c_{t}=\bar{f}\left( k_{t},k_{t+1}\right)
\end{eqnarray}

We can then write the Bellman equation for the discretized problem as%
\begin{equation}
V\left( k_{t}\right) =\max_{k_{t+1}\in K}u\left( k_{t},\bar{f}\left(
k_{t},k_{t+1}\right) \right) +\beta V\left( k_{t+1}\right) .
\end{equation}

As you see, this is a Bellman for a constrained problem, i.e., the control
variable is constrained relative to the case when $k_{t+1}$ is continuous.
Two things should be noted; First, the Bellman equation is the true Bellman
equation of the constrained problem and previous results hold, in particular
the contraction mapping theorems apply. Second, how important the constraint
implied by the discretization depends on how fine the grid is. Many (few)
elements of $K$ with small (large) distances between them, imply that the
constraint is weak (severe).

Denoting an arbitrary initial value function by $V_{0}\left( k_{t}\right) ,$
being $n$ numbers, we update this value function according to%
\begin{equation}
V_{1}\left( k_{t}\right) =\max_{k_{t+1}\in K}u\left( k_{t},\bar{f}\left(
k_{t},k_{t+1}\right) \right) +\beta V_{0}\left( k_{t+1}\right)
\end{equation}%
giving $V_{1}\left( k_{t}\right) ,$ being a new set of $n$ numbers. We then
iterate on the Bellman equation%
\begin{equation}
V_{s+1}\left( k_{t}\right) =\max_{k_{t+1}\in K}u\left( k_{t},\bar{f}\left(
k_{t},k_{t+1}\right) \right) +\beta V_{s}\left( k_{t+1}\right)
\end{equation}%
until $V_{s+1}\left( k_{t}\right) \approx V_{s}\left( k_{t}\right) .$ For
each of the $n$ values of $k_{t}$, we check $u\left( k_{t},\bar{f}\left(
k_{t},k_{t+1}\right) \right) +\beta V\left( k_{t+1}\right) $ for all $n$
values of $k_{t+1}$ and choose the $k_{t+1}$ that gives the highest value,
giving $V_{s+1}\left( k_{t}\right) $. Therefore, each iteration requires $%
n^{2}$ evaluations when the state variable is unidimensional.\footnote{%
When the state variable is higher of dimensionality, this method quickly
becomes computationally too burdensome.}

Lets consider a simple example. 
\begin{eqnarray}
&&\max_{\left\{ k_{t+1},c_{t}\right\} _{t}^{\infty }}\sum_{t=0}^{\infty
}\beta ^{t}\ln \left( c_{t}\right) \\
\text{s.t. }k_{t+1} &=&f\left( k_{t},c_{t}\right) \equiv k_{t}^{\alpha
}+\left( 1-\delta \right) k_{t}-c_{t} \\
k_{t} &\geq &0\forall t \\
k_{0} &=&\underline{k}
\end{eqnarray}

First, we solve for 
\begin{equation}
c_{t}=\bar{f}\left( k_{t},k_{t+1}\right) =k_{t}^{\alpha }+\left( 1-\delta
\right) k_{t}-k_{t+1}.
\end{equation}

Then, we note that $k_{t}\in \left[ 0,\delta ^{\frac{1}{\alpha -1}}\right]
\implies $ $k_{t+1}\in \left[ 0,\delta ^{\frac{1}{\alpha -1}}\right] $,
implying that value function is bounded. If also $\beta <1,$ the Bellman
equation%
\begin{equation}
V\left( k_{t}\right) =\max_{c_{t}}\ln \left( k_{t}^{\alpha }+\left( 1-\delta
\right) k_{t}-k_{t+1}\right) +\beta V\left( k_{t+1}\right)
\end{equation}%
is a contraction mapping.

Let us parametrize, setting $\beta =0.9$, $\delta =.2$ and $\alpha
=1/2\implies \delta ^{\frac{1}{\alpha -1}}=25$ and discretize the state
space by requiring 
\begin{equation}
k_{t}\in \left[ 5,10,15,20,25\right] \equiv K\text{~}\forall t.
\end{equation}%
Now set an initial value function, for example 
\begin{equation*}
V_{0}\left( k\right) =\ln k~\forall k.
\end{equation*}

This is then updated in the following way. For each possible $k_{t},k_{t+1}$
we calculate the left hand side of the Bellman equation, and solve the
maximization problem.\ So, for $k_{t}=5,$ and all $k_{t+1}\in K$ we have%
\begin{eqnarray*}
\ln \left( 5^{\alpha }+\left( 1-\delta \right) 5-5\right) +.9\ln 5 &=&1.66 \\
\ln \left( 5^{\alpha }+\left( 1-\delta \right) 5-10\right) +.9\ln 10
&=&-\infty \\
\ln \left( 5^{\alpha }+\left( 1-\delta \right) 5-15\right) +.9\ln 15
&=&-\infty \\
\ln \left( 5^{\alpha }+\left( 1-\delta \right) 5-20\right) +.9\ln 20
&=&-\infty \\
\ln \left( 5^{\alpha }+\left( 1-\delta \right) 5-25\right) +.9\ln 25
&=&-\infty
\end{eqnarray*}

Implying that the updated value function for $k_{t}=5$ is 
\begin{equation*}
V_{1}\left( 5\right) =1.66.
\end{equation*}%
For $k_{t}=10,$%
\begin{eqnarray*}
\ln \left( 10^{\alpha }+\left( 1-\delta \right) 10-5\right) +.9\ln 5 &=&3.27
\\
\ln \left( 10^{\alpha }+\left( 1-\delta \right) 10-10\right) +.9\ln 10
&=&2.22 \\
\ln \left( 10^{\alpha }+\left( 1-\delta \right) 10-15\right) +.9\ln 15
&=&-\infty \\
\ln \left( 10^{\alpha }+\left( 1-\delta \right) 10-20\right) +.9\ln 20
&=&-\infty \\
\ln \left( 10^{\alpha }+\left( 1-\delta \right) 10-25\right) +.9\ln 25
&=&-\infty
\end{eqnarray*}%
implying 
\begin{equation*}
V_{1}\left( 10\right) =3.27
\end{equation*}

In the same way, for $k_{t}=15$ 
\begin{eqnarray*}
\ln \left( 15^{\alpha }+\left( 1-\delta \right) 15-5\right) +.9\ln 5 &=&3.83
\\
\ln \left( 15^{\alpha }+\left( 1-\delta \right) 15-10\right) +.9\ln 10
&=&3.84 \\
\ln \left( 15^{\alpha }+\left( 1-\delta \right) 15-15\right) +.9\ln 15
&=&2.30 \\
\ln \left( 15^{\alpha }+\left( 1-\delta \right) 15-20\right) +.9\ln 20
&=&-\infty \\
\ln \left( 15^{\alpha }+\left( 1-\delta \right) 15-5\right) +.9\ln 25
&=&-\infty
\end{eqnarray*}
\begin{equation*}
V_{1}\left( 15\right) =3.84
\end{equation*}

Doing this also for $k_{t}=20$ and $k_{t}=25$ completes the first iteration.
Then, we repeat the iterations until we think the process has converged
sufficiently,

\subsection{Band Matrix Methods for differential equations}

Assume we want to solve the differential equation 
\begin{equation}
y^{\prime \prime }\left( t\right) +ay^{\prime }\left( t\right) +by\left(
t\right) =g\left( t\right)
\end{equation}%
over some interval, with initial conditions given. If we want to solve this
numerically, we first have to get rid of the abstract infinitely small
differences called differentials. We approximate these with finite size
forward differences such that 
\begin{eqnarray}
y^{\prime }\left( t\right) &\approx &\frac{y\left( t+\Delta t\right)
-y\left( t\right) }{\Delta t}, \\
y^{\prime }\left( t\right) &\approx &\frac{\frac{y\left( t+\Delta t\right)
-y\left( t\right) }{\Delta t}-\frac{y\left( t\right) -y\left( t-\Delta
t\right) }{\Delta t}}{\Delta t} \\
&=&\frac{y\left( t+\Delta t\right) -2y\left( t\right) +y\left( t-\Delta
t\right) }{\Delta t^{2}}
\end{eqnarray}

Using this we can solve the equation for a finite set of values in the
following way. Say we want to solve the equation in the interval $t\in
\lbrack p,q]$ and we know $y^{\prime \prime }\left( 0\right) =c_{1}$ and $%
y^{\prime }\left( 0\right) =c_{2}.$ We divide the interval for $t$ into $n$
equal parts and use the following notation 
\begin{eqnarray}
t_{k} &=&p+\frac{k\left( q-p\right) }{n}, \\
t_{k}-t_{k-1} &=&\frac{q-p}{n}\equiv \Delta t, \\
y\left( t_{k}\right) &\equiv &y_{k.}
\end{eqnarray}

This gives us the following equations 
\begin{eqnarray}
\frac{\left( y_{-1}-2y_{0}+y_{1}\right) }{\Delta t^{2}} &=&c_{1} \\
\frac{\left( -y_{0}+y_{1}\right) }{\Delta t} &=&c_{2} \\
\frac{\left( y_{-1}-2y_{0}+y_{1}\right) }{\Delta t^{2}}+a\frac{-y_{0}+y_{1}}{%
\Delta t}+by_{0} &=&g\left( t_{0}\right) \\
\frac{\left( y_{0}-2y_{1}+y_{2}\right) }{\Delta t^{2}}+a\frac{-y_{1}+y_{2}}{%
\Delta t}+by_{1} &=&g\left( t_{1}\right) \\
&&. \\
&&. \\
\frac{\left( y_{n-1}-2y_{n}+y_{n+1}\right) }{\Delta t^{2}}+a\frac{%
-y_{n}+y_{n+1}}{\Delta t}+by_{n} &=&g\left( t_{n}\right) .
\end{eqnarray}

This provides $n+3$ linear equations for the $n+3$ unknown $y$. Writing this
as a system we have 
\begin{eqnarray}
\mathbf{A}\left[ 
\begin{array}{c}
y_{-1} \\ 
y_{0} \\ 
y_{1} \\ 
. \\ 
y_{n} \\ 
y_{n+1}%
\end{array}%
\right] &=&\left[ 
\begin{array}{c}
c_{1} \\ 
c_{2} \\ 
g\left( t_{0}\right) \\ 
g\left( t_{1}\right) \\ 
. \\ 
g\left( t_{n}\right)%
\end{array}%
\right] \\
\left[ 
\begin{array}{c}
y_{-1} \\ 
y_{0} \\ 
y_{1} \\ 
. \\ 
y_{n} \\ 
y_{n+1}%
\end{array}%
\right] &=&\mathbf{A}^{-1}\left[ 
\begin{array}{c}
c_{1} \\ 
c_{2} \\ 
g\left( t_{0}\right) \\ 
g\left( t_{1}\right) \\ 
. \\ 
g\left( t_{n}\right)%
\end{array}%
\right] ,
\end{eqnarray}%
with (setting $n=3$) 
\begin{gather}
\mathbf{A\equiv } \\
_{\left[ 
\begin{array}{cccccc}
_{\Delta t^{-2}} & _{-2\Delta t^{-2}} & _{\Delta t^{-2}} & _{0} & _{0} & _{0}
\\ 
_{0} & _{-\Delta t^{-1}} & _{\Delta t^{-1}} & _{0} & _{0} & _{0} \\ 
_{\Delta t^{-2}} & _{-2\Delta t^{-2}-a\Delta t^{-1}+b} & _{\Delta
t^{-2}+a\Delta t^{-1}} & _{0} & _{0} & _{0} \\ 
_{0} & _{\Delta t^{-2}} & _{-2\Delta t^{-2}-a\Delta t^{-1}+b} & _{\Delta
t^{-2}+a\Delta t^{-1}} & _{0} & _{0} \\ 
_{0} & _{0} & _{\Delta t^{-2}} & _{-2\Delta t^{-2}-a\Delta t^{-1}+b} & 
_{\Delta t^{-2}+a\Delta t^{-1}} & _{0} \\ 
_{0} & _{0} & _{.} & _{\Delta t^{-2}} & _{-2\Delta t^{-2}-a\Delta t^{-1}+b}
& _{\Delta t^{-2}+a\Delta t^{-1}}%
\end{array}%
\right] }  \notag
\end{gather}

To get any accuracy, we should of course set $n$ much larger than $3.$ As we
see, the matrix $\mathbf{A}$ contains many zeros, with a \emph{band }of
non-zeros around the diagonal. Due to this feature, it is easy for the
computer to invert it also if $n$ is in the order of hundreds.

\subsection{Newton-Raphson}

Suppose we are looking for an optimum of the real valued function $f\left( 
\mathbf{x}\right) ,\mathbf{x\in R}^{n}$ where $f$ is twice differentiable. A
standard way to do this numerically is to apply the Newton-Raphson
algorithm. If the optimum is interior, it satisfies the necessary first
order conditions that the gradient is zero, i.e., at the optimum, denoted $%
\mathbf{x}^{\ast }$, 
\begin{equation}
Df\left( \mathbf{x}\right) =0.
\end{equation}

Now apply a first order linear approximation to the gradient from some
initial point $\mathbf{x}_{0}$ 
\begin{equation}
0=Df\left( \mathbf{x}^{\ast }\right) \approx Df\left( \mathbf{x}_{0}\right)
+D^{2}f\left( \mathbf{x}_{0}\right) \left( \mathbf{x}^{\ast }-\mathbf{x}%
_{0}\right) ,
\end{equation}
where $D^{2}f\left( \mathbf{x}_{0}\right) $ is the \emph{Hessian} matrix of
second derivatives of $f.$

Provided the Hessian is invertible, we can get an approximation to $\mathbf{x%
}^{\ast },$%
\begin{equation}
\mathbf{x}^{\ast }\approx \mathbf{x}_{0}-\left( D^{2}f\left( \mathbf{x}%
_{0}\right) \right) ^{-1}Df\left( \mathbf{x}_{0}\right) .
\end{equation}

From this we can construct a search algorithm that under some circumstances
makes better and better approximations 
\begin{equation}
\mathbf{x}_{s+1}\approx \mathbf{x}_{s}-\left( D^{2}f\left( \mathbf{x}%
_{s}\right) \right) ^{-1}Df\left( \mathbf{x}_{s}\right) .
\end{equation}

If we don't have analytic solutions to the gradient and Hessian we can use
numerical approximations, for example the forward difference method; For a
small number $\varepsilon $, we have, 
\begin{eqnarray}
\frac{\partial f\left( \mathbf{x}\right) }{\partial x_{1}} &=&\frac{f\left( 
\mathbf{x}+\left[ 
\begin{array}{c}
\varepsilon \\ 
\mathbf{0}%
\end{array}%
\right] \right) -f\left( \mathbf{x}\right) }{\varepsilon } \\
\frac{\partial f^{2}\left( \mathbf{x}\right) }{\partial x_{1}^{2}} &=&\frac{%
f\left( \mathbf{x}+\left[ 
\begin{array}{c}
\varepsilon \\ 
\mathbf{0}%
\end{array}%
\right] \right) -2f\left( \mathbf{x}\right) +f\left( \mathbf{x-}\left[ 
\begin{array}{c}
\varepsilon \\ 
\mathbf{0}%
\end{array}%
\right] \right) }{\varepsilon ^{2}}, \\
\frac{\partial f^{2}\left( \mathbf{x}\right) }{\partial x_{2}\partial x_{1}}
&=&\frac{\frac{\partial f\left( \mathbf{x}+\left[ 
\begin{array}{c}
0 \\ 
\varepsilon \\ 
\mathbf{0}%
\end{array}%
\right] \right) }{\partial x_{1}}-\frac{\partial f\left( \mathbf{x}\right) }{%
\partial x_{1}}}{\varepsilon } \\
&=&\frac{f\left( \mathbf{x}+\left[ 
\begin{array}{c}
\varepsilon \\ 
\varepsilon \\ 
\mathbf{0}%
\end{array}%
\right] \right) -f\left( \left[ 
\begin{array}{c}
0 \\ 
\varepsilon \\ 
\mathbf{0}%
\end{array}%
\right] \right) -f\left( \mathbf{x}+\left[ 
\begin{array}{c}
\varepsilon \\ 
0 \\ 
\mathbf{0}%
\end{array}%
\right] \right) +f\left( \mathbf{x}\right) }{\varepsilon ^{2}}.
\end{eqnarray}

One should be very careful with this method since it can only find local
optima in the case when $f$ is not globally concave. In well-behaved
problems, it is however, easily programmed and fairly quick. \newpage

\section{Note on constant hazard models}

%TCIMACRO{\TeXButton{TeX field}{\setcounter{equation}{0}} }%
%BeginExpansion
\setcounter{equation}{0}
%EndExpansion
Consider as an example a pool of unemployed, with measure (size) at time $t$
given by $x(t).$ Suppose also that there is an outflow of people from the
unemployment pool (we are for now disregarding the inflow of new
unemployed). The outflow is determined by an assumption that there is a
constant probability per unit of time, denoted $h$ to get hired. Using this
we can derive a law-of-motion for $x\left( t\right) $. Over a small
(infinitesimal) interval of time $dt$, we have 
\begin{eqnarray}
x\left( t+dt\right) &=&x\left( t\right) \left( 1-hdt\right) \\
x\left( t+dt\right) -x\left( t\right) &=&-x\left( t\right) hdt \\
\lim_{dt\rightarrow 0}\frac{x\left( t+dt\right) -x\left( t\right) }{dt}
&\equiv &\dot{x}\left( t\right) =-hx\left( t\right) .
\end{eqnarray}

This is a simple differential equation with a solution 
\begin{equation}
x\left( t\right) =x\left( 0\right) e^{-ht}.  \label{poolsize}
\end{equation}

Now, consider an individual who is unemployed at time $0$ and the random
variable $s$ denote the time she will stay unemployed. Let us now derive the
probability density function $f\left( s\right) $ and let $F\left( s\right) $
denote the cumulative distribution function, i.e., $F\left( s\right) $ is
the probability the unemployment spell is no longer than $s.$ Clearly, 
\begin{equation}
F\left( s\right) =\int_{0}^{s}f\left( t\right) dt.
\end{equation}

From (\ref{poolsize}) we know that at time $s$, a share 
\begin{equation}
\frac{x\left( t\right) }{x\left( 0\right) }=e^{-hs},
\end{equation}

remains unemployed and since hiring is completely random, 
\begin{equation}
1-\frac{x\left( t\right) }{x\left( 0\right) }=1-e^{-hs}=F\left( s\right) ,
\end{equation}

and consequently 
\begin{equation}
f\left( s\right) =he^{-hs}.
\end{equation}

Let us also verify that the probability of finding a job per unit of time,
conditional on not having found it stays constant. This probability is 
\begin{equation}
\frac{f\left( t\right) }{1-F\left( t\right) }=\frac{he^{-ht}}{e^{-ht}}=h.
\end{equation}

We can now compute the average spell length as 
\begin{equation}
\int_{0}^{\infty }sf\left( s\right) ds=\int_{0}^{\infty }she^{-hs}ds.
\end{equation}

Using the formula for integration by parts, 
\begin{eqnarray}
\int_{0}^{\infty }she^{-hs}ds &=&\left[ -se^{-hs}\right] _{0}^{\infty
}-\int_{0}^{\infty }-e^{-hs}ds \\
&=&0-0-\left[ \frac{e^{-hs}}{h}\right] _{0}^{\infty }=\frac{1}{h}.
\end{eqnarray}

Similarly, we can compute the median length $m$, i.e., solving $F(m)=1/2$
from 
\begin{eqnarray}
1/2 &=&1-e^{-hm} \\
&\rightarrow &e^{-hm}=1/2 \\
-hm &=&-\ln 2 \\
m &=&\frac{\ln 2}{h}\approx \frac{.69}{h}
\end{eqnarray}

This is sometimes called the \emph{rule of 69; }expressing $h$ in percent
per unit of time. The half-life is found by dividing $69$ by $h.$ For
example, if the probability of finding job is 5\% per week. It takes $%
69/5\approx 14$ weeks before half the pool of unemployed have found jobs and
the average unemployment spell is $1/0.05=20$ weeks.

Another example; it is often found that differences in GDP between countries
(after controlling for differences is savings and schooling) is closing by $%
3\%$ per year. Then, half the difference is left after $69/3=23$. The rule
of 69 also works when something is growing. If a bank account yields 4\%
return per year, it takes $69/4\approx 17$ years for it to double.

\section{T-mappings\label{Tmap}}

Instead of using time as an argument of the value function, let's use time
subscripts. We can then write the Bellman equation as 
\begin{equation*}
V_{t}\left( k_{t}\right) =\left\{ 
\begin{array}{c}
\max_{c_{t}}\left( u\left( k_{t},c_{t}\right) +\beta V_{t+1}\left(
k_{t+1}\right) \right) \\ 
s.t.k_{t+1}=f\left( k_{t},c_{t}\right)%
\end{array}%
\right.
\end{equation*}

Now, let us define an operator that maps next period value functions$%
,V_{t+1}:K\rightarrow R,$ (where $K$ is the state space, i.e., the set of
possible values for the state variable), into functions that provides the
current value associated with all $k_{t}$ in the state space. Thus, the
operator, which we will call,\thinspace $T,$ maps elements of the space of
value functions, call that space $C,$ back into the same space, i.e., $%
T:C\rightarrow C$.\footnote{%
Later on, we will make assumptions such that we can further restrict the
space of possible valute functions.}

Formally, we define the $T$ mapping as 
\begin{equation*}
TV_{t+1}:C\rightarrow C\equiv \left\{ 
\begin{array}{c}
\max_{c_{t}}\left( u\left( k_{t},c_{t}\right) +\beta V_{t+1}\left(
k_{t+1}\right) \right) \\ 
s.t.k_{t+1}=f\left( k_{t},c_{t}\right)%
\end{array}%
\right. .
\end{equation*}%
When we want to indicate that the mapped function, $TV_{t+1},$ is a function
of, e.g., $k_{t},$we append $\left( k_{t}\right) ;$%
\begin{equation*}
TV_{t+1}=TV_{t+1}\left( k_{t}\right) .
\end{equation*}

Note that while $V_{t+1}$ is a function of $k_{t+1}$, $TV_{t+1}$ is a
function of $k_{t}.$

Let us take an example. Suppose $u\left( k_{t},c_{t}\right) =\ln c_{t}$ and $%
f\left( k_{t},c_{t}\right) =k_{t}-c_{t}$ $\forall t.$ Let's us now see what
the $T$ operator does. Take a particular element in $C,$ for example $\ln
k:R^{+}\rightarrow R.$ So here the state space $K=R^{+}.$ Now, 
\begin{equation*}
T\ln k_{t}=\left\{ 
\begin{array}{c}
\max_{c_{t}}\left( \ln c_{t}+\beta \left( \ln k_{t+1}\right) \right) \\ 
s.t.k_{t+1}=k_{t}-c_{t}%
\end{array}%
\right.
\end{equation*}

The first order condition is 
\begin{equation*}
\frac{1}{c_{t}}=\frac{\beta }{k_{t}-c_{t}}\rightarrow c_{t}=\frac{k_{t}}{%
1+\beta }
\end{equation*}%
thus, 
\begin{eqnarray*}
T\ln k_{t} &=&\ln \frac{k_{t}}{1+\beta }+\beta \left( \ln \left( k_{t}-\frac{%
k_{t}}{1+\beta }\right) \right) \\
&=&\left( \left( 1+\beta \right) \ln k_{t}+\beta \ln \beta -\left( 1+\beta
\right) \ln \left( 1+\beta \right) \right)
\end{eqnarray*}
which is another function $R^{+}\rightarrow R.$ So, we see that the $T$ maps
functions into (potentially) other functions, concluding the example.

Using the $T$ operator, the Bellman equation in period $t$ can be written 
\begin{equation*}
V_{t}=TV_{t+1},
\end{equation*}%
or, equivalently, 
\begin{equation*}
V_{t}\left( k_{t}\right) =TV_{t+1}\left( k_{t}\right) .
\end{equation*}

Furthemore, in the next period, the next period's Bellman equation is 
\begin{equation*}
V_{t+1}=TV_{t+2}.
\end{equation*}

Thus, 
\begin{equation*}
V_{t}=TV_{t+1}=T^{2}V_{t+2}.
\end{equation*}

The meaning of $T^{2}V_{t+2}$ in words is; give me (I am $T^{2}$) a value
function that applies in period $t+2$ (you don't need to say anything about
what $k_{t+2}$ is going to be), and I give you a function that tells you the
value in period $t$ associated with all values of $k_{t}$ in the state
space. Formally;%
\begin{eqnarray*}
T^{2}V_{t+2} &=&\left\{ 
\begin{array}{c}
\max_{c_{t}}\left( u\left( k_{t},c_{t}\right) +\beta \left( TV_{t+2}\left(
k_{t+1}\right) )\right) \right) \\ 
s.t.k_{t+1}=f\left( k_{t},c_{t}\right)%
\end{array}%
\right. \\
&=&\left\{ 
\begin{array}{c}
\max_{c_{t}}\left( u\left( k_{t},c_{t}\right) +\beta \left(
\max_{c_{t+1}}\left( u\left( k_{t+1},c_{t+1}\right) +\beta V_{t+2}\left(
k_{t+2}\right) \right) \right) \right) \\ 
s.t.k_{t+1}=f\left( k_{t},c_{t}\right) ,k_{t+2}=f\left(
k_{t+1},c_{t+1}\right)%
\end{array}%
\right. .
\end{eqnarray*}%
Now, suppose in the time autonomuous case that we can find a limiting
function 
\begin{equation*}
V\left( k_{t}\right) \equiv \lim_{s\rightarrow \infty }T^{s}V_{t+s}\left(
k_{t}\right)
\end{equation*}%
then, as shown in the lecture notes, this function satisfies the Bellman
equation 
\begin{equation*}
V\left( k_{t}\right) =TV\left( k_{t}\right)
\end{equation*}%
i.e., it is a fixed point of the $T$ operator. This means that in the the
space $C$, $V$ is an elements such that $T$ maps back onto the same element.
If $T$ is a contraction mapping on $C$, this element exists and is unique.

\end{document}
