Saturday, December 24, 2011

Compuational economics lecture 23

Posted by Unknown at 11:18 PM


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)

 

0 comments:

Post a Comment

Subscribe to: Post Comments (Atom)