Our basic dynamic model

An arbitrary infinite horizon problem can be represented using a Bellman equation $$V(S_t) = \max_{q_t} u(q_t) + \beta\,V(S_{t+1}(q_t))$$

Our basic dynamic model

An arbitrary infinite horizon problem can be represented using a Bellman equation $$V(S_t) = \max_{q_t} u(q_t) + \beta\,V(S_{t+1}(q_t))$$

We can start from any arbitrary initial state vector and simulate the optimal policy paths in deterministic or stochastic settings

Our basic dynamic model

An arbitrary infinite horizon problem can be represented using a Bellman equation $$V(S_t) = \max_{q_t} u(q_t) + \beta\,V(S_{t+1}(q_t))$$

We can start from any arbitrary initial state vector and simulate the optimal policy paths in deterministic or stochastic settings

We know that this problem has a fixed point (solution) $V^*$ from previous lectures, but how do we approximate $V^*$?

Our basic dynamic model

An arbitrary infinite horizon problem can be represented using a Bellman equation $$V(S_t) = \max_{q_t} u(q_t) + \beta\,V(S_{t+1}(q_t))$$

We can start from any arbitrary initial state vector and simulate the optimal policy paths in deterministic or stochastic settings

We know that this problem has a fixed point (solution) $V^*$ from previous lectures, but how do we approximate $V^*$?

One way to do this is via projection

Our basic dynamic model

An arbitrary infinite horizon problem can be represented using a Bellman equation $$V(S_t) = \max_{q_t} u(q_t) + \beta\,V(S_{t+1}(q_t))$$

We can start from any arbitrary initial state vector and simulate the optimal policy paths in deterministic or stochastic settings

We know that this problem has a fixed point (solution) $V^*$ from previous lectures, but how do we approximate $V^*$?

One way to do this is via projection

Or perturbation, or discretization, or some analytic approaches for specific models

Projection methods

Main idea: build some function $\hat{V}$ indexed by coefficients that approximately solves the Bellman

Projection methods

Main idea: build some function $\hat{V}$ indexed by coefficients that approximately solves the Bellman

What do I mean by approximately?

Projection methods

Main idea: build some function $\hat{V}$ indexed by coefficients that approximately solves the Bellman

What do I mean by approximately?

The coefficients of $\hat{V}$ are selected to minimize some residual function
that tells us how far away our solution is from solving the Bellman

Projection methods

Main idea: build some function $\hat{V}$ indexed by coefficients that approximately solves the Bellman

What do I mean by approximately?

The coefficients of $\hat{V}$ are selected to minimize some residual function
that tells us how far away our solution is from solving the Bellman

Rearrange the Bellman equation and define a new functional $H$ that maps the problem into a more general framework $$H(V) = V(S_t) - \max_{q_t} u(q_t) + \beta\,V(S_{t+1}(q_t))$$

Projection methods

How do we do this?

Projection methods

We then define a residual $$R(\mathbf{S}|\theta) = H(V^j(\mathbf{S}|\theta))$$ and select the coefficient values to minimize the residual, given some measure of distance

Projection methods

We then define a residual $$R(\mathbf{S}|\theta) = H(V^j(\mathbf{S}|\theta))$$ and select the coefficient values to minimize the residual, given some measure of distance

This step of selecting the coefficients is called projecting $H$ against our basis

Projection methods

At first this may be a pretty confusing concept

Projection methods

At first this may be a pretty confusing concept

Simple example to get intuition?

Projection methods

At first this may be a pretty confusing concept

Simple example to get intuition?

Ordinary least squares linear regression

Projection methods

We don't know the true functional form, but we can approximate it using the first two monomials on $X$: $1$ and $X$ $$E[Y|X] \approx \theta_0 + \theta_1 X$$

Projection methods

We don't know the true functional form, but we can approximate it using the first two monomials on $X$: $1$ and $X$ $$E[Y|X] \approx \theta_0 + \theta_1 X$$

These are the first two elements of the monomial basis

Projection methods

We don't know the true functional form, but we can approximate it using the first two monomials on $X$: $1$ and $X$ $$E[Y|X] \approx \theta_0 + \theta_1 X$$

These are the first two elements of the monomial basis

One residual function is then: $R(Y,X|\theta_0,\theta_1) = abs(Y-\theta_0-\theta_1 X)$

Projection methods

We don't know the true functional form, but we can approximate it using the first two monomials on $X$: $1$ and $X$ $$E[Y|X] \approx \theta_0 + \theta_1 X$$

These are the first two elements of the monomial basis

One residual function is then: $R(Y,X|\theta_0,\theta_1) = abs(Y-\theta_0-\theta_1 X)$

In econometrics we often minimize the square of this residual, which is just OLS

Projection methods

We don't know the true functional form, but we can approximate it using the first two monomials on $X$: $1$ and $X$ $$E[Y|X] \approx \theta_0 + \theta_1 X$$

These are the first two elements of the monomial basis

One residual function is then: $R(Y,X|\theta_0,\theta_1) = abs(Y-\theta_0-\theta_1 X)$

In econometrics we often minimize the square of this residual, which is just OLS

So OLS is within the class of projection methods

Projection classes are defined by metrics

Projection methods are separated into several broad classes by the type of residual we're trying to shrink to zero

Projection classes are defined by metrics

Projection methods are separated into several broad classes by the type of residual we're trying to shrink to zero

We need to select some metric function $\rho$, that determines how we project

Example residuals given different projections

Example: The figure shows two different residuals on some capital domain of $[0,\bar{k}]$

The residual based on the coefficient vector $\theta_1$ is large for small values of capital but near-zero everywhere else

Example residuals given different projections

Example: The figure shows two different residuals on some capital domain of $[0,\bar{k}]$

The residual based on the coefficient vector $\theta_1$ is large for small values of capital but near-zero everywhere else

The residual based on $\theta_2$ has medium values just about everywhere

We want to use weighted residuals

We move from the plain residual to $\rho$ because we want to set a weighted residual equal to zero (more general)

We want to use weighted residuals

We move from the plain residual to $\rho$ because we want to set a weighted residual equal to zero (more general)

Suppose we have some weight functions $\phi_i:\Omega \rightarrow \mathbb{R}$ that map from our state space to the real line

We want to use weighted residuals

We move from the plain residual to $\rho$ because we want to set a weighted residual equal to zero (more general)

Suppose we have some weight functions $\phi_i:\Omega \rightarrow \mathbb{R}$ that map from our state space to the real line

The one-dimensional metric is defined as $$\begin{equation*} \rho(R\cdot|\theta,0) = \begin{cases} 0 & \text{if } \int_\Omega\phi_i(\mathbf{S})R(\cdot|\theta)d\mathbf{S} = 0, i=1,...,j+1\\ 1 &\text{otherwise} \end{cases} \end{equation*}$$

We want to use weighted residuals

We can then change our problem to simply solving a system of integrals ensuring the metric is zero $$\int_\Omega\phi_i(\mathbf{S})R(\cdot|\theta)d\mathbf{S} = 0, i=1,...,j+1.$$

We want to use weighted residuals

We can then change our problem to simply solving a system of integrals ensuring the metric is zero $$\int_\Omega\phi_i(\mathbf{S})R(\cdot|\theta)d\mathbf{S} = 0, i=1,...,j+1.$$

We can solve this using standard rootfinding techniques

We want to use weighted residuals

We can then change our problem to simply solving a system of integrals ensuring the metric is zero $$\int_\Omega\phi_i(\mathbf{S})R(\cdot|\theta)d\mathbf{S} = 0, i=1,...,j+1.$$

We can solve this using standard rootfinding techniques

Big remaining question: how do we choose our $j+1$ weight functions?

Least squares projection

Suppose we selected the weight function to be $$\phi_i(\mathbf{S}) = \frac{\partial R(\mathbf{S}|\theta)}{\partial \theta_{i-1}}$$

Least squares projection

Recall the objective of least squares is $$\min_{\mathbf{\theta}} \int R^2(\cdot|\theta)d\mathbf{S}$$

Least squares projection

Recall the objective of least squares is $$\min_{\mathbf{\theta}} \int R^2(\cdot|\theta)d\mathbf{S}$$ The FOC for a minimum is $$\int \frac{\partial R(\mathbf{S}|\theta)}{\partial \theta_{i-1}}R(\cdot|\theta)d\mathbf{S} = 0, i = 1,...,j+1$$

Least squares projection

Recall the objective of least squares is $$\min_{\mathbf{\theta}} \int R^2(\cdot|\theta)d\mathbf{S}$$ The FOC for a minimum is $$\int \frac{\partial R(\mathbf{S}|\theta)}{\partial \theta_{i-1}}R(\cdot|\theta)d\mathbf{S} = 0, i = 1,...,j+1$$

So the first order condition sets the weighted average residual to zero
where the weights are determined by the partial derivatives of the residual

Collocation

The most commonly used weight function gives us a methodology called collocation

Collocation

The most commonly used weight function gives us a methodology called collocation

Here our weight function is $$\phi_i(\mathbf{S}) = \delta(\mathbf{S}-\mathbf{S}_i)$$ Where $\delta$ is the Dirac delta function and $\mathbf{S}_i$ are the j+1 collocation points or collocation nodes selected by the researcher

Collocation

The most commonly used weight function gives us a methodology called collocation

Here our weight function is $$\phi_i(\mathbf{S}) = \delta(\mathbf{S}-\mathbf{S}_i)$$ Where $\delta$ is the Dirac delta function and $\mathbf{S}_i$ are the j+1 collocation points or collocation nodes selected by the researcher

The Dirac delta function is zero at all $\mathbf{S}$ except at $\mathbf{S} = \mathbf{S}_i$

Collocation

Before we even select our coefficients, this means that the residual can only be non-zero at a finite set of points $\mathbf{S_i}$

Collocation

Before we even select our coefficients, this means that the residual can only be non-zero at a finite set of points $\mathbf{S_i}$

So the solution to our problem must set the residual to zero at these collocation points

Rough idea of how we proceed

Collocation methods aim to approximate the value function $V(S_t)$ with some linear combination of known, and usually orthogonal $(\int f(x)g(x)w(x)=0)$, basis functions

Rough idea of how we proceed

Collocation methods aim to approximate the value function $V(S_t)$ with some linear combination of known, and usually orthogonal $(\int f(x)g(x)w(x)=0)$, basis functions

One example of a possible class of basis functions is the monomial basis: $x, x^2, x^3,...$

Rough idea of how we proceed

Recall we solve for coefficients $\theta$ by setting the residual to be zero at all of our collocation points

Rough idea of how we proceed

Recall we solve for coefficients $\theta$ by setting the residual to be zero at all of our collocation points

But the residual is a function of $V$, and thus will be a function of $\theta$, a loop seemingly impossible to escape

Rough idea of how we proceed

In any given iteration, we:

Rough idea of how we proceed

In any given iteration, we:

1. Begin with a vector of coefficients on the basis functions (e.g. a random guess)

Rough idea of how we proceed

In any given iteration, we:

1. Begin with a vector of coefficients on the basis functions (e.g. a random guess)

2. Use a linear combination of the basis functions as an approximation to the value function on the right hand side of the Bellman

Rough idea of how we proceed

In any given iteration, we:

1. Begin with a vector of coefficients on the basis functions (e.g. a random guess)

2. Use a linear combination of the basis functions as an approximation to the value function on the right hand side of the Bellman

3. Solve the Bellman with the approximated value function in it, at our set of collocation points

Rough idea of how we proceed

In any given iteration, we:

1. Begin with a vector of coefficients on the basis functions (e.g. a random guess)

2. Use a linear combination of the basis functions as an approximation to the value function on the right hand side of the Bellman

3. Solve the Bellman with the approximated value function in it, at our set of collocation points

4. Recover a set maximized continuation values at these points in our state space conditional on the previous value function approximant we used

Still have some questions

By the contraction mapping theorem this is guaranteed to converge from any arbitrary initial set of coefficients! (conditional on satisfying the assumptions and there being no numerical error..)

Still have some questions

By the contraction mapping theorem this is guaranteed to converge from any arbitrary initial set of coefficients! (conditional on satisfying the assumptions and there being no numerical error..)

This still leaves several questions unanswered

Still have some questions

By the contraction mapping theorem this is guaranteed to converge from any arbitrary initial set of coefficients! (conditional on satisfying the assumptions and there being no numerical error..)

This still leaves several questions unanswered

• Where in the state space do we place the collocation points?

Still have some questions

By the contraction mapping theorem this is guaranteed to converge from any arbitrary initial set of coefficients! (conditional on satisfying the assumptions and there being no numerical error..)

This still leaves several questions unanswered

• Where in the state space do we place the collocation points?

• What basis functions do we use?

Still have some questions

By the contraction mapping theorem this is guaranteed to converge from any arbitrary initial set of coefficients! (conditional on satisfying the assumptions and there being no numerical error..)

This still leaves several questions unanswered

• Where in the state space do we place the collocation points?

• What basis functions do we use?

Do your answers to the previous two questions really matter?

Still have some questions

Schemes exist to generate high quality approximations, and to obtain the approximation at low computational cost

Still have some questions

Schemes exist to generate high quality approximations, and to obtain the approximation at low computational cost

Once we have recovered an adequate approximation to the value function, we have effectively solved the entire problem!

Still have some questions

Schemes exist to generate high quality approximations, and to obtain the approximation at low computational cost

Once we have recovered an adequate approximation to the value function, we have effectively solved the entire problem!

We know how the policymaker's expected discounted stream of future payoffs changes as we move through the state space

Still have some questions

Schemes exist to generate high quality approximations, and to obtain the approximation at low computational cost

Once we have recovered an adequate approximation to the value function, we have effectively solved the entire problem!

We know how the policymaker's expected discounted stream of future payoffs changes as we move through the state space

We can solve the Bellman at some initial state and know that solving the Bellman will yield us optimal policies

Interpolation

What points in our state space do we select to be collocation points?

Interpolation

What points in our state space do we select to be collocation points?

We often have continuous states in economics (capital, temperature, oil reserves, technology shocks, etc.), so we must have some way to reduce the infinity of points in our state space into something more manageable

Interpolation

Using our knowledge of how the value function behaves at the limited set of points on our grid, we can interpolate our value function approximant at all points off the grid points, but within the domain of our grid

Basis functions

Let $V$ be the value function we wish to approximate with some $\hat{V}$

Basis functions

Let $V$ be the value function we wish to approximate with some $\hat{V}$

$\hat{V}$ is constructed as a linear combination of $n$ linearly independent (i.e. orthogonal) basis functions $$\hat{V}(x)= \sum_{j=1}^n c_j \psi_j(x)$$

Basis functions

The number of basis functions we select, $n$, is the degree of interpolation

Basis functions

The number of basis functions we select, $n$, is the degree of interpolation

In order to recover $n$ coefficients, we need at least $n$ equations that must be satisfied at a solution to the problem

Basis functions

The number of basis functions we select, $n$, is the degree of interpolation

In order to recover $n$ coefficients, we need at least $n$ equations that must be satisfied at a solution to the problem

If we have precisely $n$ equations, we are just solving a simple system of linear equations: we have a perfectly identified system and are solving a collocation problem

Basis functions

Solve a system of equations, linear in $c_j$ that equates the function approximant at the grid points to the recovered values $$\Psi \mathbf{c} = \mathbf{y}$$ where $\Psi$ is the matrix of basis functions, $c$ is a vector of coefficients, and $y$ is a vector of the recovered values

Basis functions

Solve a system of equations, linear in $c_j$ that equates the function approximant at the grid points to the recovered values $$\Psi \mathbf{c} = \mathbf{y}$$ where $\Psi$ is the matrix of basis functions, $c$ is a vector of coefficients, and $y$ is a vector of the recovered values

We can recover $c$ by left dividing by $\Psi$ which yields $$\mathbf{c} = \Psi^{-1}\mathbf{y}$$

Basis functions

Solve a system of equations, linear in $c_j$ that equates the function approximant at the grid points to the recovered values $$\Psi \mathbf{c} = \mathbf{y}$$ where $\Psi$ is the matrix of basis functions, $c$ is a vector of coefficients, and $y$ is a vector of the recovered values

We can recover $c$ by left dividing by $\Psi$ which yields $$\mathbf{c} = \Psi^{-1}\mathbf{y}$$

If we have more equations, or grid points, than coefficients, then we can just use OLS to solve for the coefficients by minimizing the sum of squared errors

(Pseudo-)spectral methods

Spectral methods apply all of our basis functions to the entire domain of our grid: they are global

(Pseudo-)spectral methods

Spectral methods apply all of our basis functions to the entire domain of our grid: they are global

When using spectral methods we virtually always use polynomials

(Pseudo-)spectral methods

Spectral methods apply all of our basis functions to the entire domain of our grid: they are global

When using spectral methods we virtually always use polynomials

Why?

(Pseudo-)spectral methods

Spectral methods apply all of our basis functions to the entire domain of our grid: they are global

When using spectral methods we virtually always use polynomials

Why?

The Stone-Weierstrass Theorem which states (for one dimension)

Suppose $f$ is a continuous real-valued function defined on the interval $[a,b]$.
For every $\epsilon > 0, \,\, \exists$ a polynomial $p(x)$ such that for all $x \in [a,b]$ we have $||f(x) - p(x)||_{sup} \leq \epsilon$

(Pseudo-)spectral methods

For any continuous function $f(x)$, we can approximate it arbitrarily well with some polynomial $p(x)$, as long as $f(x)$ is continuous

(Pseudo-)spectral methods

For any continuous function $f(x)$, we can approximate it arbitrarily well with some polynomial $p(x)$, as long as $f(x)$ is continuous

This means the function can even have kinks

(Pseudo-)spectral methods

For any continuous function $f(x)$, we can approximate it arbitrarily well with some polynomial $p(x)$, as long as $f(x)$ is continuous

This means the function can even have kinks

Unfortunately we do not have infinite computational power so solving kinked problems with spectral methods is not advised

Basis choice

What would be your first choice of basis?

Basis choice

What would be your first choice of basis?

Logical choice: the monomial basis: $1, x, x^2,...$

Basis choice

What would be your first choice of basis?

Logical choice: the monomial basis: $1, x, x^2,...$

It is simple, and SW tells us that we can uniformly approximate any continuous function on a closed interval using them

Basis choice

What would be your first choice of basis?

Logical choice: the monomial basis: $1, x, x^2,...$

It is simple, and SW tells us that we can uniformly approximate any continuous function on a closed interval using them

Practice

code up a function project_monomial(f, n, lb, ub) that takes in some function f, degree of approximation n, lower bound lb and upper bound ub, and constructs a monomial approximation on an evenly spaced grid via collocation

Approximating sin(x)

Now we need to construct a function f_approx(coefficients, plot_points) that takes in the coefficients vector, and an arbitrary vector of points to evaluate the approximant at (for plotting)

Approximating sin(x)

Now we need to construct a function f_approx(coefficients, plot_points) that takes in the coefficients vector, and an arbitrary vector of points to evaluate the approximant at (for plotting)

function f_approx(coefficients, points)
    n = length(coefficients) - 1
    basis_functions = [coefficients[degree + 1] * points.^degree for degree = 0:n] # evaluate basis functions
    basis_matrix = hcat(basis_functions...)                                        # transform into matrix
    function_values = sum(basis_matrix, dims = 2)                                  # sum up into function value
    return function_values
end;

Cool!

We just wrote some code that exploits Stone-Weierstrauss and allows us to (potentially) approximate any continuous function arbitrarily well as n goes to infinity

Cool!

We just wrote some code that exploits Stone-Weierstrauss and allows us to (potentially) approximate any continuous function arbitrarily well as n goes to infinity

To approximate any function we'd need to feed in some basis function g(x, n) as opposed to hard-coding it like I did in the previous slides

Cool!

We just wrote some code that exploits Stone-Weierstrauss and allows us to (potentially) approximate any continuous function arbitrarily well as n goes to infinity

To approximate any function we'd need to feed in some basis function g(x, n) as opposed to hard-coding it like I did in the previous slides

Turns out we never use the monomial basis though

Cool!

Turns out we never use the monomial basis though

Cool!

Turns out we never use the monomial basis though

Why?

Monomials are not good

What's the deal?

Monomials are not good

What's the deal?

The matrix of monomials, $\Phi$, is often ill-conditioned, especially as the degree of the monomials increases

Monomials are not good

What's the deal?

The matrix of monomials, $\Phi$, is often ill-conditioned, especially as the degree of the monomials increases

The first 6 monomials can induce a condition number of $10^{10}$, a substantial loss of precision

Monomials are not good

runge(x) = 1 ./ (1 .+ 25x.^2);
coefficients_10 = project_monomial(runge, 10, -1, 1);
points = rand(10);
n = length(coefficients_10) - 1;
basis_functions = [coefficients_10[degree + 1] * points.^degree for degree = 0:n];
basis_matrix = hcat(basis_functions...);
println("Condition number: $(cond(basis_matrix))")

## Condition number: 1.7447981077501208e20

Monomials are not good

runge(x) = 1 ./ (1 .+ 25x.^2);
coefficients_10 = project_monomial(runge, 10, -1, 1);
points = rand(10);
n = length(coefficients_10) - 1;
basis_functions = [coefficients_10[degree + 1] * points.^degree for degree = 0:n];
basis_matrix = hcat(basis_functions...);
println("Condition number: $(cond(basis_matrix))")

## Condition number: 1.7447981077501208e20

The Chebyshev basis

Most frequently used is the Chebyshev basis

The Chebyshev basis

Most frequently used is the Chebyshev basis

It has nice approximation properties:

1. They are easy to compute

The Chebyshev basis

Most frequently used is the Chebyshev basis

It has nice approximation properties:

1. They are easy to compute

2. They are orthogonal

The Chebyshev basis

Most frequently used is the Chebyshev basis

It has nice approximation properties:

1. They are easy to compute

2. They are orthogonal

3. They are bounded between $[-1,1]$

The Chebyshev basis

Most frequently used is the Chebyshev basis

It has nice approximation properties:

1. They are easy to compute

2. They are orthogonal

3. They are bounded between $[-1,1]$

Chebyshev polynomials are often selected because they minimize the oscillations that occur when approximating functions like Runge's function

The Chebyshev basis

Chebyshev polynomials are defined by a recurrence relation, $$\begin{gather} T_0(x) = 1 \\ T_1(x) = x \\ T_{n+1} = 2xT_n(x) - T_{n-1}(x) \end{gather}$$ and are defined on the domain $[-1,1]$

The Chebyshev basis

Chebyshev polynomials are defined by a recurrence relation, $$\begin{gather} T_0(x) = 1 \\ T_1(x) = x \\ T_{n+1} = 2xT_n(x) - T_{n-1}(x) \end{gather}$$ and are defined on the domain $[-1,1]$

In practice this is easy to expand to any interval $[a,b]$

The Chebyshev basis

Chebyshev polynomials are defined by a recurrence relation, $$\begin{gather} T_0(x) = 1 \\ T_1(x) = x \\ T_{n+1} = 2xT_n(x) - T_{n-1}(x) \end{gather}$$

Write two functions: cheb_polys(n, x) and monomials(n, x) with a degree of approximation n and vector of points x, that return the values of the approximant at x

The Chebyshev basis

Chebyshev polynomials are defined by a recurrence relation, $$\begin{gather} T_0(x) = 1 \\ T_1(x) = x \\ T_{n+1} = 2xT_n(x) - T_{n-1}(x) \end{gather}$$

Write two functions: cheb_polys(n, x) and monomials(n, x) with a degree of approximation n and vector of points x, that return the values of the approximant at x

Write a plotting function plot_function(basis_function, x, n) that takes an arbitrary basis function basis_function and plots all basis functions up to degree n

Chebyshev polynomials rule

Chebyshev polynomials span the space

Chebyshev polynomials rule

Chebyshev polynomials span the space

Monomials clump together

Chebyshev polynomials rule

Chebyshev polynomials span the space

Monomials clump together

Chebyshev polynomials are nice for approximation because they are orthogonal and they form a basis with respect to the weight function $\phi(x) = 1/\sqrt{1-x^2}$

Chebyshev polynomials rule

Chebyshev polynomials span the space

Monomials clump together

Chebyshev polynomials are nice for approximation because they are orthogonal and they form a basis with respect to the weight function $\phi(x) = 1/\sqrt{1-x^2}$

They span the polynomial vector space

Chebyshev polynomials rule

Chebyshev polynomials span the space

Monomials clump together

Chebyshev polynomials are nice for approximation because they are orthogonal and they form a basis with respect to the weight function $\phi(x) = 1/\sqrt{1-x^2}$

They span the polynomial vector space

This means that you can form any polynomial of degree equal to less than the Chebyshev polynomial you are using

Chebyshev zeros and alternative rep

Also note that the Chebyshev polynomial of order $n$ has $n$ zeros given by $$x_k = cos\left(\frac{2k-1}{2n}\pi \right), \,\,\, k=1,...,n$$ which tend to cluster quadratically towards the edges of the domain

Two important theorems

There are two important theorems to know about Chebyshev polynomials

Two important theorems

There are two important theorems to know about Chebyshev polynomials

Chebyshev interpolation theorem: If $f(x) \in \mathbb{C}[a,b]$, if $\{\psi_i(x), i=0,...\}$ is a system of polynomials (where $\psi_i(x)$ is of exact degree i) orthogonal with respect to $\phi(x)$ on $[a,b]$ and if $p_j = \sum_{i=0}^j \theta_i \psi_i(x)$ interpolates $f(x)$ in the zeros of $\psi_{n+1}(x)$, then: $$\lim_{j\rightarrow\infty} \left(|| f-p_j||_2 \right)^2 = \lim_{j\rightarrow\infty}\int_a^b \phi(x) \left(f(x) - p_j \right)^2 dx = 0$$

Two important theorems

If we have an approximation set of basis functions that are exact at the roots of the $n^{th}$ order polynomials, then as $n$ goes to infinity the approximation error becomes arbitrarily small and converges at a quadratic rate

Two important theorems

If we have an approximation set of basis functions that are exact at the roots of the $n^{th}$ order polynomials, then as $n$ goes to infinity the approximation error becomes arbitrarily small and converges at a quadratic rate

This holds for any type of polynomial, but if they are Chebyshev then convergence is uniform

Two important theorems

Chebyshev truncation theorem: The error in approximating $f$ is bounded by the sum of all the neglected coefficients

Two important theorems

Chebyshev truncation theorem: The error in approximating $f$ is bounded by the sum of all the neglected coefficients

Since Chebyshev polynomials satisfy Stone-Weierstrauss, an infinite sequence of them can perfectly approximate any continuous function

Two important theorems

Chebyshev truncation theorem: The error in approximating $f$ is bounded by the sum of all the neglected coefficients

Since Chebyshev polynomials satisfy Stone-Weierstrauss, an infinite sequence of them can perfectly approximate any continuous function

Since Chebyshev polynomials are bounded between $[-1,1]$, the sum of the omitted terms is bounded by the sum of the magnitude of the coefficients

Two important theorems

We often also have that Chebyshev approximations geometrically converge which give us the following convenient property: $$|f(x) - f^j(x|\theta)| \sim O(\theta_j)$$ The truncation error by stopping at polynomial $j$ is of the same order as the magnitude of the coefficient $\theta_j$ on the last polynomial

Boyd's moral principle

Chebyshev polynomials are the most widely used basis

Boyd's moral principle

Chebyshev polynomials are the most widely used basis

This is not purely theoretical but also from practical experience

Grid point selection

We construct the basis function approximant by evaluating the function on a predefined grid in the domain of $V$

Grid point selection

We can write this problem more compactly as $$\Phi c = y \tag{interpolation equation}$$ where

• $y$ is the column vector of $V(x_i)$
• $c$ is the column vector of coefficients $c_j$
• $\Phi$ is an $n\times n$ matrix of the $n$ basis functions evaluated at the $n$ points

Grid point selection

We can write this problem more compactly as $$\Phi c = y \tag{interpolation equation}$$ where

• $y$ is the column vector of $V(x_i)$
• $c$ is the column vector of coefficients $c_j$
• $\Phi$ is an $n\times n$ matrix of the $n$ basis functions evaluated at the $n$ points

If we recover a set of values at our interpolation nodes, $V^*(x_i)$, we can then simply invert $\Phi$ and right multiply it by $y$ to recover our coefficients

Chebyshev strikes again

A good selection of points are called Chebyshev nodes

Chebyshev strikes again

A good selection of points are called Chebyshev nodes

These are simply the roots of the Chebyshev polynomials above on the domain $[-1,1]$

Chebyshev strikes again

A good selection of points are called Chebyshev nodes

These are simply the roots of the Chebyshev polynomials above on the domain $[-1,1]$

They are given by $$x_k = cos\left(\frac{2k-1}{2n}\pi\right),\,\, k = 1,...,n$$ for some Chebyshev polynomial of degree $n$

Chebyshev strikes again

A good selection of points are called Chebyshev nodes

These are simply the roots of the Chebyshev polynomials above on the domain $[-1,1]$

They are given by $$x_k = cos\left(\frac{2k-1}{2n}\pi\right),\,\, k = 1,...,n$$ for some Chebyshev polynomial of degree $n$

Mathematically, these also help reduce error in our approximation

Chebyshev node locations

The nodes are heavily focused near the end points, why is that?

Chebyshev node locations

The nodes are heavily focused near the end points, why is that?

Imagine areas of our approximant near the center of our domain but not at a node

Chebyshev node locations

The nodes are heavily focused near the end points, why is that?

Imagine areas of our approximant near the center of our domain but not at a node

These areas benefit from having multiple nodes on both the left and right

Chebyshev node locations

The nodes are heavily focused near the end points, why is that?

Imagine areas of our approximant near the center of our domain but not at a node

These areas benefit from having multiple nodes on both the left and right

This provides more information for these off-node areas and allows them to be better approximated because we know whats happening nearby in several different directions

Chebyshev node locations

The nodes are heavily focused near the end points, why is that?

Imagine areas of our approximant near the center of our domain but not at a node

These areas benefit from having multiple nodes on both the left and right

This provides more information for these off-node areas and allows them to be better approximated because we know whats happening nearby in several different directions

If we moved to an area closer to the edge of the domain, there may only one node to the left or right of it providing information on what the value of our approximant should be

Discrete states

How do we handle a discrete state $S_d$ when trying to approximate $V$?

Multi-dimensional approximation

Thus far we have displayed the Chebyshev basis in only one dimension

Multi-dimensional approximation

Thus far we have displayed the Chebyshev basis in only one dimension

We approximate functions of some arbitrary dimension $N$ by taking the tensor of vectors of the one-dimensional Chebyshev polynomials

Multi-dimensional approximation

Thus far we have displayed the Chebyshev basis in only one dimension

We approximate functions of some arbitrary dimension $N$ by taking the tensor of vectors of the one-dimensional Chebyshev polynomials

Construct a vector of polynomials $[\phi_{1,1}, \, \phi_{1,2}, \, \phi_{1,3}]$ for dimensions 1

Multi-dimensional approximation

Thus far we have displayed the Chebyshev basis in only one dimension

We approximate functions of some arbitrary dimension $N$ by taking the tensor of vectors of the one-dimensional Chebyshev polynomials

Construct a vector of polynomials $[\phi_{1,1}, \, \phi_{1,2}, \, \phi_{1,3}]$ for dimensions 1

Construct a vector of polynomials $[\phi_{2,1}, \, \phi_{2,2}, \, \phi_{2,3}]$ for dimension 2

Multi-dimensional approximation

Thus far we have displayed the Chebyshev basis in only one dimension

We approximate functions of some arbitrary dimension $N$ by taking the tensor of vectors of the one-dimensional Chebyshev polynomials

Construct a vector of polynomials $[\phi_{1,1}, \, \phi_{1,2}, \, \phi_{1,3}]$ for dimensions 1

Construct a vector of polynomials $[\phi_{2,1}, \, \phi_{2,2}, \, \phi_{2,3}]$ for dimension 2

The tensor is just the product of every possibly polynomial pair which results in: $[\phi_{1,1}\phi_{2,1} ,\, \phi_{1,1}\phi_{2,2}, \, \phi_{1,1}\phi_{2,3}, \, \phi_{1,2}\phi_{2,1}, \, \phi_{1,2}\phi_{2,2}, \, \phi_{1,2}\phi_{2,3}, \, \phi_{1,3}\phi_{2,1}, \, \phi_{1,3}\phi_{2,2}, \, \phi_{1,3}\phi_{2,3}]$

Multi-dimensional approximation

The computational complexity here grows exponentially: $\text{total # points} = (\text{points per state})^{\text{# states}}$

Multi-dimensional approximation

The computational complexity here grows exponentially: $\text{total # points} = (\text{points per state})^{\text{# states}}$

Exponential complexity is costly, often called the Curse of dimensionality

Finite element methods

An alternative to pseudo-spectral methods is the class of finite element methods

Finite element methods

An alternative to pseudo-spectral methods is the class of finite element methods

Finite element methods use basis functions that are non-zero over subintervals of the domain of our grid

Finite element methods

The higher the order the polynomial we use, the higher the order of derivatives that we can preserve continuity at the knots

Finite element methods

The higher the order the polynomial we use, the higher the order of derivatives that we can preserve continuity at the knots

For example, a linear spline yields an approximant that is continuous, but its first derivatives are discontinuous step functions unless the underlying value function happened to be precisely linear

Finite element methods

As we increase the order of the spline polynomial, we have increasing numbers of coefficients we need to determine

Finite element methods

As we increase the order of the spline polynomial, we have increasing numbers of coefficients we need to determine

To determine these additional coefficients using the same number of points, we require additional conditions that must be satisfied

Finite element methods

With linear splines, each segment of our value function approximant is defined by a linear function

Finite element methods

With linear splines, each segment of our value function approximant is defined by a linear function

For each of these linear components, we need to solve for 1 coefficient and 1 intercept term

Finite element methods

With linear splines, each segment of our value function approximant is defined by a linear function

For each of these linear components, we need to solve for 1 coefficient and 1 intercept term

Each end of a linear segment must equal the function value at the knots

Finite element methods

With linear splines, each segment of our value function approximant is defined by a linear function

For each of these linear components, we need to solve for 1 coefficient and 1 intercept term

Each end of a linear segment must equal the function value at the knots

We have two conditions and two unknowns for each segment: this is a simple set of linear equations that we can solve

Cubic splines

Suppose we wish to approximate using a cubic spline on $N+1$ knots

Cubic splines

Suppose we wish to approximate using a cubic spline on $N+1$ knots

We need $N$ cubic polynomials when entails $4N$ coefficients to determine

Cubic splines

Suppose we wish to approximate using a cubic spline on $N+1$ knots

We need $N$ cubic polynomials when entails $4N$ coefficients to determine

We can obtain $3(N-1)$ equations by ensuring that the approximant is continuous at all interior knots, and its first and second derivatives are continuous at all interior knots $[3\times(N+1-1-1)]$

Cubic splines

Ensuring the the approximant equals the function's value at all of the nodes adds another $N+1$ equations

Cubic splines

Ensuring the the approximant equals the function's value at all of the nodes adds another $N+1$ equations

We therefore have a total of $4N-2$ equations for $4N$ coefficients

Cubic splines

Ensuring the the approximant equals the function's value at all of the nodes adds another $N+1$ equations

We therefore have a total of $4N-2$ equations for $4N$ coefficients

We need two more conditions to solve the problem

Splines can potentially handle lots of curvature

If the derivative is of interest for optimization, or to recover some variable of economic meaning, then we may need to have these derivatives preserved well at the knots

Splines can potentially handle lots of curvature

If the derivative is of interest for optimization, or to recover some variable of economic meaning, then we may need to have these derivatives preserved well at the knots

One large benefit of splines is that they can handle kinks or areas of high curvature

Splines can potentially handle lots of curvature

If the derivative is of interest for optimization, or to recover some variable of economic meaning, then we may need to have these derivatives preserved well at the knots

One large benefit of splines is that they can handle kinks or areas of high curvature

How?

Splines can potentially handle lots of curvature

If the derivative is of interest for optimization, or to recover some variable of economic meaning, then we may need to have these derivatives preserved well at the knots

One large benefit of splines is that they can handle kinks or areas of high curvature

How?

By having the modeler place many knots in a concentrated region

Splines can potentially handle lots of curvature

If the derivative is of interest for optimization, or to recover some variable of economic meaning, then we may need to have these derivatives preserved well at the knots

One large benefit of splines is that they can handle kinks or areas of high curvature

How?

By having the modeler place many knots in a concentrated region

If the knots are stacked on top of one another, this actually allows for a kink to be explicitly represented in the value function approximant, but the economist must know precisely where the kink is