Computational Economics

Go to top, next, previous, or johnstachurski.net

Function Approximation

A major topic when we need to work with functions

Important for dynamic programming

The Problem with Finite Memory

In numerical computing we often work with arrays

Take an array,
apply some operation to get a new array,
apply another operation to get yet another array, etc.

At each stage we can store the modified array, because it only requires finite memory

Example:

Returning to the topic of Markov chains (see this lecture)

Let p be a Markov matrix
- represented as a 2D NumPy array
Let q be a distribution (probability mass function)
- represented as a one dimensional NumPy array

If the current state X[t] has distribution q, then X[t+1] has distribution np.dot(q, p)

See EDTC, section 4.2.2 for details

So if q is the distribution of the initial state X[0], then the time distribution of X[t] can be computed via

for j in range(t):
    q = np.dot(q, p)   # The final value of q is the distribution of X[t]

At each iteration we only need to store the finite array q, so there's no problem here

But now suppose we want to work with functions, rather than arrays

Consider this function:

def f(x):
    return np.sin(x) - 2 * np.cos(x)  # Assume numpy imported as np

This function f is easily stored on the computer

Defined in terms of finite number of elementary functions, constants, etc.

However, suppose we want to

Take a mathematical function,
apply some operation to get a new function,
apply another operation to get yet another function, etc.

In practice this might be problematic, because functions are inherently infinite

Can't store the newly created functions at each step

Example:

In the next lecture we will study value function iteration

Begin with a function w
Create a new function w' according to

where rho, f, U and phi are known functions/constants

The number of possible y is infinite, so we can't evaluate at every y in finite time

Instead we approximate w' as follows

Evaluate w' on a finite grid of points, and store these values
Build an approximation to w' using this data

More generally, to approximate a function we need

A finite amount of data giving information about the function
A rule which tells us how to reconstruct the function (approximately) from the data

This is naturally done using OOP, since we are combining data and behavior (reconstruction of the function)

We will concentrate on simple approximation schemes suitable for dynamic programming

Step Functions

One of the simplest approximation schemes is step functions

Useful because

Sometimes simplicity is a good thing
Interact well with dynamic programming (see EDTC, section 6.2.2)

A one-dimensional step function s is defined by

an increasing array X (the grid), and
an array of values Y with len(Y) = len(X)

The function s is the (broken) blue line in this figure

A simplistic implementation of s is as follows

def s(x):
    if x < X[0]:      # x < first element of grid
        y = 0
    elif x >= X[-1]:  # x >= last element of grid
        y = Y[-1]
    else:
        for i in range(len(X)-1):
            if X[i] <= x < X[i+1]:
                y = Y[i]
    return y

Sometimes we want to evaluate the expectation of s(Z) where

Z is a random variable
The cumulative distribution function of Z is given by F

For example, the expectation of the simple function in the figure is

Y[0] * (F(X[1]) - F(X[0])) + Y[1] * (F(X[2]) - F(X[1])) + Y[2] * (1 - F(X[2]))

Problem:

Design a class StepFun such that an instance represents a particular step function

The instance variables are X and Y
The class implements a __call__ method that evaluates s(x) for any x
The class implements an expectation method that
- takes as an argument a cumulative distribution function F
- and returns the expectation of s(Z) when F is the distribution of X

Try to be efficient (i.e., avoid loops)

Hint: If X is an increasing array, then X.searchsorted(x, side='right') - 1 returns the largest i such that X[i] <= x

Solution:

import numpy as np

class StepFun:

    def __init__(self, X, Y):
        """
        Parameters: 

            * X is an increasing array of length n
            * Y is an array of length n

        Implements a step function s, where 
        
            s(x) = sum_{i=0}^{n-1} Y[i] 1{X[i] <= x < X[i+1]}

        where X[n] := infty
        """
        self.X, self.Y = X, Y

    def __call__(self, x):
        "Evaluate the step function at x."
        if x < self.X[0]:
            return 0.0
        i = self.X.searchsorted(x, side='right') - 1
        return self.Y[i]

    def expectation(self, F):
        "Computes expectation of s(Y) when F is the cdf of Y."
        probs = np.append(F(self.X), 1)
        return np.dot(self.Y, probs[1:] - probs[:-1])

Linear Interpolation

Now we turn to another approximation scheme: piecewise linear interpolation

Linear interpolation is implemented by the NumPy/SciPy interp function

Here's an example of usage

import numpy as np
import matplotlib.pyplot as plt
xvals = np.linspace(-2, 2, 100)
yvals = np.sin(xvals)
plt.plot(xvals, yvals, 'k-')
eval_points = [-2, -1, 0, 1, 2]
plt.plot(eval_points, np.interp(eval_points, xvals, yvals), 'r-')
plt.show()

Generates this graph

When we do dynamic programming, we will use a simple class which wraps the function interp

## Filename: lininterp.py
## Author: John Stachurski

from scipy import interp

class LinInterp:
    "Provides linear interpolation in one dimension."

    def __init__(self, X, Y):
        """Parameters: X and Y are sequences or arrays
        containing the (x,y) interpolation points."""
        self.X, self.Y = X, Y

    def __call__(self, z):
        return interp(z, self.X, self.Y)

The purpose of the class is to store the interpolation data locally (i.e., as instance data)

Lectures on Computational Economics

John Stachurski

Function Approximation

The Problem with Finite Memory

Step Functions

Linear Interpolation