In squipbar/cheby: Computes Chebychev approximations to 1- and 2-dimensional functions

Introduction

The cheby package provides easy-to-use Chebychev interpolation of an arbitrary R function.

Chebychev interpolation approximates an arbitrary function on a closed interval using an the Chebychev polynomials as a basis. In general, polynomial approximation is helpful when a function $f:[a,b]\rightarrow \mathbb R$ is expensive to compute at all values in the domain. So instead, one approximates $f(\cdot)$ is by $\hat f(\cdot)$, where:

$$\hat f(x) = \sum_{i=1}^n a_i \phi_i(x) $$

Where $\phi_i(x)$ are polynomials. In practice, this involves computing the function at a limited number of grid points, and then fitting a polynomial function through the function values at these points. This fitting is typically done by linear regression of the points on the fitting polynomials. This minimizes the square deviations at the grid points, but imposes no restrictions on the shape of the approximation.

If we want the resulting approximation to have particular properties (such as monotonicity or concavity), then we will need to use other methods. This package provides functionality to fit one-dimensional functions both via conventional Chebychev approximation, and by imposing constraints on the shape of the approximation.

The algorithms in this package draw heavily on Judd (1998).

One-dimensional Chebychev approximation

The Theory

Algorithm 6.2 of Judd (1998) outlines the steps required to produce the one dimensional Chebychev interpolation of order $n$ on $m$ nodes (with $n<m$) of the function $f:[a,b]\rightarrow \mathbb R$.

1. Compute the interpolation nodes on $[-1,1]$ via: $$z_k = -\cos \left( \frac{2k-1}{2m}\pi \right)$$
2. Map the nodes to the interval $[a,b]$ $$x_k = \left( \frac{b-a}{2} \right)(z_k+1) + a$$
3. Compute $f(\cdot)$ at the nodes $$y_k = f(x_k), \qquad k=1,\ldots,m$$
4. Regress the $y_k$ on the Chebychev polynomials $T_i(z)$ for $i=0,\ldots,n$ $$a_i = \frac{\sum_{k=1}^m y_k T_i(z_k) }{\sum_{k=1}^m T_i(z_k)^2}$$
5. The Chebychev approximation is then given by: $$\hat f(x) = \sum_{i=1}^n a_i T_i\left(\frac{2(x-1)}{b-a}-1\right)$$

Theorem 6.7.3 of Judd (1998) shows that this approximation is pointwise-convergent for any $C^k$ function.

Basic approximation

The basic command for Chebychev interpolation is d1.poly. For example, to compute the approximation of the natural logarithm:

suppressMessages(library(cheby))
base <- d1.poly( fn=log, range=c(.01,4), iOrder=10, iPts=11 )

This computes a 6th-order approximation to log on the interval $[0.01,4]$ (recall that log is not continuous at zero) using 7 grid points (the minimum possible). The output is a function that we can use straight away. We can check the approximation at a few points:

c( base(1), log(1) )
c( base(2), log(2) )

So the approximation is good to 2 decimal places for these values of $x$. Better still, we can see the approximation visually:

plot( log, xlim=c(.01, 4), type='l', lty=2 ) 
plot( base, xlim=c(.01, 4), col='red', add=TRUE, lwd=2 )
legend( 'bottomright', c('log', 'Order 10 polynomial approximation'), 
        col=c('black', 'red'), lty=2:1, bty='n', lwd=1:2 )

So the approximation is very good for most of the approximation range, although given the discontinuity at zero it struggles for small $x$. And it is fast:

library(microbenchmark)
microbenchmark(base <- d1.poly( fn=log, range=c(.01,4), iOrder=6, iPts=7 ))

However, the approximation is not strictly concave - it oscillates slightly around the approximated function. Note how the dashed log function shows up alternately above and below the approximation. In some numerical calculations, preserving shape properties of the function can be important. Later we will discuss how value function will often fail if the continuation value function is not concave, as then there is no unique local interior maximum. This is not a problem which can be eliminated by simply increasing the order of the approximation, either. In fact this makes the problem worse. A higher order approximation improves the fit for small $x$ but comes at the cost of more noticeable oscillations at large $x$.

Extracting more information about the approximation

To extract more information about the approximation, set the details flag to TRUE.

base.deets <- d1.poly( fn=log, range=c(.01,4), iOrder=10, iPts=11, details=TRUE )

This returns a five-element list with elements:

• fn The polynomial approximation function. So `base.deets$fn(1)$ will equal $base(1)$
• poly The polynomial representation of the Chebychev approximation. This element is a polynomial from the polynom package, so can be manipulated easily (see the documentation of that package for more details). It is also defined only on the $[-1,1]$ interval, so needs to be translated to the approximation interval to match the output of the fn element.
• fn.deriv The exact derivative of the polynomial approximation function. Computed by differentiating the polynomial representation and translating to the approximation interval. This is very useful for optimization problems where the derivatives are required. Plotting this function (not shown here) will generate upward-sloping parts of the function, confirming that the approximation is not globally concave. We can check also this calculation by comparing to the finite-difference calculation:

1e06 * ( base(1+1e-06) - base(1) )
base.deets$fn.deriv(1)

• fn.deriv.poly The polynomial representation of the derivative. The analogue of poly, but for deriv.
• residuals The residuals of the approximation. All zero if iPts is one more than iOrder.

Pre-computing function values

The "standard" usage of d1.poly above integrates function calculation and approximation. This is convenient - it simply allows one to pass a function to the approximation and get a function back. But it is not very flexible - it requires a self-contained function to be evaluated separately at each of the points in the grid. In many applications this is not efficient. For example, if the target function to be evaluated contains a computationally demanding common component, it may be quicker to calculate that on its own and then include that value in the evaluation of $f(\cdot)$ at each of the grid points.

Here is an example of how to use this syntax

grid.pts <- d1.grid( c( .01, 4), 11 )
    # The grid of Chebychev nodes
log.vals <- sapply( grid.pts, log )
    # The target function evaluated at the approximation nodes
base.grid <- d1.poly( fn.vals=log.vals, grid=grid.pts, range=c(.01, 4), iOrder=10, iPts=11, details=TRUE )

The inclusion of the argument grid here is superfluous, as without this d1.poly assumes that the function values are evaluated on a grid of Chebychev nodes, which is exactly what d1.grid produces here. If you wish to evaluate the function on a different grid of nodes, you should use this argument to specify that grid.

Passing function parameters

One can pass function parameters using the argument fn.opts. However, this requires that the target function have all its parameters enter through a single argument named opts (using a list to pass multiple parameters where required).

For example, to compute a fifth-order approximation to the function $A k ^ \alpha$ where $A$ and $\alpha$ are parameters:

target <- function( k, opts ) return( opts$A  * k ^ opts$alpha ) 
    # The target function
apx.1 <- d1.poly( target, c(.001,2), 5, 6, fn.opts=list(A=1, alpha=.7) )
apx.2 <- d1.poly( target, c(.001,2), 5, 6, fn.opts=list(A=1.1, alpha=.2) )
    # The approximations
plot( function(x) target(x, opts=list(A=1, alpha=.7) ), xlim=c(.001,2), lty=2, ylab='' )
plot( function(x) target(x, opts=list(A=1.1, alpha=.2) ), xlim=c(.001,2), lty=2, add=TRUE )
plot( apx.1, xlim=c(.001,2), col='red', lwd=2, add=TRUE )
plot( apx.2, xlim=c(.001,2), col='red', lwd=2, add=TRUE )

Approximating mulit-variable functions

Functions of more variables can be approximated using dn.poly. This works much the same as d1.poly1, but currently only works for functions of two variables.

One-dimensional shape-preserving Chebychev approximation

The theory

[tbc - from Judd]

Usage

To approximate a univariate function using a shape-preserving polynomial, we can use sp1.poly. This function takes the same inputs as d1.poly, but now also requires details on the shape-preserving restrictions on the polynomial approximation. These details are provided through the arguments n.shape and sign.deriv. The former is a vector establishing the number of Chebychev nodes at which shape restrictions should be imposed. The second is a vector of positive and negative unit values determining the sign of the restriction. For example, to restrict the approximation to be increasing at 3 points and concave at 7, one would use set n.shape=c(3,7) and sign.deriv=c(1,-1).

Here's a complete example:

base.sp <- sp1.poly( fn=log, range=c(.01,4), iOrder=10, iPts=11, n.shape=c(3,21), sign.deriv=c(1,-1) )
plot( log, xlim=c(.01,4), lty=2 )
plot( base.sp, xlim=c(.01,4), add=TRUE, lwd=2, col='red' )
plot( base, xlim=c(.01,4), add=TRUE, lwd=2, col='blue' )
legend( 'bottomright', c('log', 'Shape-preserving order \n10 polynomial approximation', 'Order 10 polynomial approximation'), 
        col=c('black', 'red', 'blue'), lty=c(2,1,1), bty='n', lwd=c(1,2,2) )

Here, we impose concavity at 21 points in the approximation interval and monotonicity at only 3 (as the standard Chebychev approximation delivers near-monotonicity already). The plot above show that this generates concavity for small values of $x$, where the standard approximation oscillates much more. The cost, of course, is a poorer fit of the approximation near zero.

An example with applications in economics

I now compute an example problem where the concavity or otherwise of the approximating function makes a difference to the resulting policy function.

The problem

I solve the neoclassical growth model with CRRA preferences and Cobb-Douglas production. The state variable is $k$, the capital-labor ratio. The value function for the problem is given by:

$$ V(k) = \max_{k'>0} \left{ \frac{c^{1-\sigma}}{1-\sigma} + \beta V(k') \right} \ \text{s.t.} \qquad k' = A k^\alpha - c + (1-\delta) k$$

This example is convenient because when $\sigma=1$ (so preferences over consumption are logarithmic) and $\delta=1$, then the model has a known following analytic solution. This means that we can check the approximations relative to the exact truth in this case.

If $\sigma=\delta=1$, then the exact solution is:

$$ V(k) = B_1 + B_2 \log k $$

Where: $$ B_1 = \frac{ \alpha\beta \log (\alpha\beta) + \log(1-\alpha\beta) }{ 1 - \beta } \qquad\qquad\qquad B_2 =\frac{\alpha}{1-\alpha\beta}$$

Value function iteration

The first approach is to use value function iteration. This computes a sequence of functions $V_1, V_2, \ldots, V_N$ via: $$ V_{n+1}(k) = T( V_n )(k) $$ Where: $$ T(V)(k) = \max_{k'>0} \left{ \log c + \beta V(k') \right} \ \text{s.t.} \qquad k' = f(k) - c + (1-\delta) k$$

This is guaranteed to converge under the infinity norm as the the operator $T$ is contraction mapping.

• Include the code for the operator $T$
• Store outer loop in a separate example file
• Distinguish between the sp and non-sp cases.

Solving the first-order conditions directly

Two-dimensional approximation

References

Judd (1998), Numerical Methods in Economics, MIT Press

squipbar/cheby documentation built on May 30, 2019, 8 a.m.