## Choosing a Solution Method

Numerical analysts have special terms that they use to classify the rates at which iterative routines converge. Specifically, a sequence of iterates x(k) is said to converge to x* at a rate of order p if there is constant C > 0 such that x(fc+1) - x* < C x(k) - x* p for sufficiently large k. In particular, the rate of convergence is said to be linear if C < 1 and p 1, superlinear if 1 < p < 2, and quadratic if p 2. The asymptotic rates of convergence of the nonlinear equation solution...

## Multivariate Affine Asset Pricing Model

As the dimension of the state process increases, the no-arbitrage PDE becomes increasingly difficult to apply in practice. There are some cases, however, for which this so-called curse of dimensionality can be avoided. The most important case is the affine asset pricing model, which has been widely applied, especially in modeling interest rate and futures price term structure and in valuing European options. Suppose that the d-dimensional risk-neutral state process S can be described by an...

## [Txrk4fTx0[additional parameters

The inputs are the name of a problem file that returns the function f, the vector of time values T and the initial conditions, xo. The fourth input is an empty matrix to make the calling syntax for rk4 compatible with the CompEcon ODE solvers. The two outputs are the vector of time values (this makes rk4 compatible with Matlab's ODE solvers) and the solution values. Unlike the suite of ODE solvers provided by Matlab, rk4 is designed to be able to compute solutions for multiple initial values....

## Plotsplotresid

Here, the COMPECON library routine nodeunif is used to generate a vector splot of 500 equally spaced states spanning the approximation interval. The residual resid is then computed and plotted at the equally-spaced states. Notice that, to compute the residual, vmax is evaluated at the residual evaluation points, not the collocation nodes. However, careful inspection of the code above reveals that vmax is designed to solve the optimization problem embedded in the Bellman equation at an arbitrary...

## Government Price Controls

To solve the government price controls of Section 8.6.3 (page 8.6.3) we again use the general solver for rational expectations models, remsolve. Recall that the equilibrium is characterized by the functional complementarity condition SEyP(x(s) + y - x(x(s) + y)) - P(s - x(s)) - c(x(s)) 0 x(s) > 0, p* > P(s x(s)), x(s) > 0 p* P(s x(s)). In the notation of the general model the state variable is (s, y), the response variable is (x, a), and the shock is y. The expectation function is with a...

## [vx valmaxvfPdelta if normvvold

And a MATLAB script that performs policy iteration for the infinite horizon model is v,x valmax v,f,P,delta pstarfstar valpol x,f,P v speye n -delta pstar fstar if norm v-vold lt tol, return, end In these implementations, maxit and tol are the maximum iterations allowed and the convergence tolerance, respectively, which have specified default values that may be overridden with the routine optset as suggested by the following examples The CompEcon Toolbox also provides two utilities for...

## K ca2S f 2c pSK2 f

Which can be solved for any initial fish stock S and industry size K . To use the rk4 solver a function returning the time derivatives for the system must be supplied. function dx fdif03 t,x,flag,beta,f,del s x 1, k x 2, temp 1 beta s. k ds 1-s . s-s. k. temp dk del s. temp.2 -2 f dx ds dk As previously mentioned, the flag variable is not used but is supplied to make rk4 compatible with the ODE solvers provided by Matlab. The solver itself is called using where x0 is a matrix of starting...

## Line Search Methods

Just as was the case with rootfinding problems, it is not always best to take a full Newton step. In fact, it may be better to either stop short or move past the Newton step. If we view the Newton step as defining a search direction, performing a one-dimensional search in that direction will generally produce improved results. In practice, it is not necessary to perform a thorough search for the best point in the Newton direction. Typically, it is sufficient to assure that successive...

## Implementation of the Collocation Method

Let us now consider the practical steps that must be taken to implement the collocation method in a computer programming environment. Below, we outline the key operations using the Matlab vector processing language, presuming access to the function approximation and numerical quadrature routines contained in the COM-PECON library. The necessary steps can be implemented in virtually any other vector processing or high-level algebraic programming language, with a level of difficulty that will...

## Constrained Optimization

The simplest constrained optimization problem involves the maximization of an objective function subject to simple bounds on the choice variable According to the Karush-Kuhn-Tucker theorem, if f is differentiable on a, b , then x is a constrained maximum for f only if it solves the complementarity problem CP f', a, b 5 Conversely, if f is concave and differentiable on a,b and x solves the complementarity problem CP f ' x , a,b , then x is a constrained maximum of f if additionally f is strictly...

## Piecewise Polynomial Splines

Piecewise polynomial splines, or simply splines for short, are a rich, flexible class of functions that may be used instead of high degree polynomials to approximate a real-valued function over a bounded interval. Generally, an order k spline consists of a series of kth order polynomial segments spliced together so as to preserve continuity of derivatives of order k 1 or less. The points at which the polynomial pieces are spliced together, v1 lt v2 lt lt vp, are called the breakpoints of the...

## Broyden S Method Inverse Jacobian Update Formula

J K ' x k - x k-1 This yields the iteration rule x x f x k - f x k-1 f x . Unlike the Newton method, the secant method requires two, rather than one starting value. The secant method is graphically illustrated in Figure 3.4. The algorithm begins with the analyst supplying two distinct guesses x 0 and x 1 for the root of f. The function f is approximated using the secant line passing through x 0 ,f x 0 and x 1 , f x 1 , whose root x 2 is accepted as an improved estimate for the root of f. The...

## Fspace fundefnbastypenaborder

Here, on input, bastype is string referencing the basis function family, which can take the values 'cheb' for Chebychev polynomial basis, 'spli' for spline basis, or 'lin' for a linear spline basis with finite difference derivatives n is the vector containing the degree of approximation along each dimension a is the vector of left endpoints of interpolation intervals in each dimension b is the vector of right Figure 6.9 Approximation of 1 25x2 1 Figure 6.9 Approximation of 1 25x2 1 endpoints of...

## End end

Here, the routine getindex is used to find the index of the following period's state. Fifth, specifies the terminal value function vterm. Upon reaching period T 1, if the energy stock is positive the animal procreates with probability 1 if the stock is zero, however, the animal procreates with probability 0 vterm ones n,1 vterm 1 0 Sixth, one packs the model structure. Here, model is a structured variable whose fields contain the reward matrix, transition probability matrix, time horizon,...

## Exercises

How many Gauss-Jacobi and Gauss-Seidel iterations are required to get answers that agree with the L-U decomposition solution to four significant digits 2.2. Use the matlab function randn to generate a random 100 by 100 matrix A and a random 100-vector b. Then use the Matlab functions tic and toc compute the time needed to solve the linear equation Ax b 1, 10, and 50 times for each of the following algorithms b x U L b , computing the L-U factors of A only once using the Matlab function lu. c x...

## Numerical Differentiation

The most natural way to approximate the derivative of a univariate real-valued function f at a point x is choose a small non-zero number h and compute a one-sided f x lim f x h f x . J w h o h Thus, at least in theory, one can compute an arbitrarily accurate approximation by taking h sufficiently small. We will take up the question of how small h should be in practice, but first we address the issue of how large an error is produced using a finite difference approximant. An error bound for the...

## Newtons Method

In practice, most nonlinear problems are solved using Newton's method or one of its variants. Newton's method is based on the principle of successive linearization. Successive linearization calls for a hard nonlinear problem to be replaced with a sequence of simpler linear problems whose solutions converge to the solution of the nonlinear problem. Newton's method is typically formulated as a rootfinding technique, but may be used to solve a fixed-point problem x g x by recasting it as the...

## Complementarity Methods

Although the complementarity problem appears quite different from the ordinary rootfinding problem, it actually can be reformulated as one. In particular, x solves the complementarity problem CP , a, b if and only if it solves the rootfinding problem x min max x ,a x ,b x 0 where min and max are applied row-wise. A formal proof of the equivalence between the complementarity problem CP f, a,b and its 'minmax' rootfinding formulation x 0 is straightforward, but requires a somewhat tedious...

## Function Iteration

Function iteration is a relatively simple technique that may be used to compute a fixed-point x g x of a function from n to n. The technique is also applicable to a rootfinding problem f x 0 by recasting it as the equivalent fixed-point problem x x f x . Function iteration begins with the analyst supplying a guess x 0 for the fixed-point of g. Subsequent iterates are generated using the simple iteration rule Since g is continuous, if the iterates converge, they converge to a fixed-point of g....

## Glossary of Matlab Terms

A list of Matlab terms used in this text. For a complete list, see Matlab documentation. chol computes the Cholesky decomposition of symmetric positive definite matrix diag returns diagonal elements of matrix as a vector or creates a diagonal matrix from a vector disp displays results on screen eps machine precision largest number that, added to 1, returns 1 eye returns order n identity matrix In eye n feval evaluates function referred to by name feval f,x,y find produces index of values...