Read this to find out how to solve the dynamic program using different methods.
⇩
Overview
Once your code has called CreateSpaces() your DP model is set up and ready to be solved. To do this, apply a solution method to the model, which is an object derived from the Method class. You can apply more than one method to your problem, although some methods only work with some kinds of models.
Every Method has a Solve() function that applies the method to the DP problem. Solution methods work at the highest level of the DP environment. That is, if your DP environment involves multiple problems (because it includes fixed and random effects), the solution method will call methods to handle each specific DP problem. (It is possible to specify that only one problem be solved, and this is used to efficiently parallelize the solution of many DP problems over a cluster.)
The most flexible way to use a solution method is to create a new object of the method and then call its Solve routine in any place of the code. This approach is necessary when using a nested algorithm in which the DP solution is embedded inside some other iteration (such as estimating parameters of the model). Some solution methods have a simpler version in which the code calls a function that creates the method for the user, applies it solve routine and then deletes the method. Both possibilities are shown here:
Example: Solve MyModel with brute force value function iteration:
Option 1 (more general)
⋮
CreateSpaces();
vmax = new ValueIteration();
vmax.Volume = LOUD; //V(θ) printed when computed
vmax -> Solve();
⋮
Loop over fixed effect variables in γ, solve each problem by calling the RandomSolve() task. Reuse the solution storage space for each fixed effect to conserve on memory.
For the current fixed group, loop over random effect variables in γ. Solve each problem by calling GSolve(). Reuse value iteration space but store choice probabilities for all random effect groups.
Given the fully specified fixed vector γ (a group), which defines an element of the group space, solve the DP problem. This is where the state space \(\Theta\) is traversed applying Bellman's equation at each point.
Inner loop of GSolve
1. Run loop over bellman iterations (or other task)
2. Update check convergence/work backwards (return to 1 until finished.)
Solution Methods are Categorized As Follows (also click on Class hierarchy at the top of any page in the DDP folder (like this page).
Within these top level categories there are some derived classes for particular algorithms.
Value Function Iteration Methods
Value function iteration applies Bellman's EMax operator, \(V(\theta)\ =\ \max U() + \delta EV(\theta') \) to solve for the value function and optimal choice probabilities.
Brute Force Iteration
ValueIteration performs brute force Bellman Equation iteration.
Example
CreateSpaces();
⋮
decl vmax = new ValueIteration();
vmax -> Solve();
Keane-Wolpin Approximation
KeaneWolpin is derived from ValueIteration and approximates the solution to the value function at a subsample of the state space.
HotzMiller: Conditional Choice Probability Methods
CCP methods use estimates of conditional choice probabilities and inverts a function of them to compute (differences in) the value function without Bellman Iteration. They rely on the CCP task which smooths data on observed choices to produce an estimate of \(P^\star(\alpha;\theta)\) that does require a solution to \(V(\theta)\).
Hotz-Miller
HotzMiller is the basic CCP method that inverts a function of the estimated CCPs to compute (differences in) the value function without Bellman Iteration.
Example
CreateSpaces();
⋮
HM = new HotzMiller(data);
Aguirregabiria and Mira
AguirregabiriaMira is derived from HotzMiller. It iterates on the Hotz-Miller inverstion to gain efficiency in estimates
Reservation values
ReservationValues is a method for solving for cut-off values in a one-dimensional continuous state variable. The distribution of the variable and the way it enters utility is unrestricted, but there are three restrictions:
The continuous vector \(\zeta\) must contain only one variable: \(\zeta = (z)\), not the typical vector of action-specific shocks.
The action vector \(\alpha\) must contain only one action: \(\alpha = (a)\). This does not mean a binary choice: a can take on any number of values. To enforce this, MyModel must be derived from the OneDimensionalChoice class.
The value of actions must satisfy the reservation property:
The model is a continuous version of the model in GetStarted.
Finite time horizon (T=10) to search over job offers
Offers are continuous: \(z \sim \Phi(x)\).
Wages not prices: \(U = (1-a)\eta + az\)
Optimal choice: a reservation wage z* at each t. That is: policy iteration not value function iteration.
Optimal Behaviour: accept offers above a reservation value z*.
z* solves a non-linear system:
$$v(1,z*)-v(0,z*)=0$$
Required Information to solve z*
4 functions of z*
Utility differences at any candidate z*: \(U(0,z*)-U(1,z*)\). This information is received by the user-provided method called Uz(z). It returns a column vector with U(0,z*)|U(1,z*)|….. (More generally when there are more N choices and N-1 cut-off values it should return a N × N-1 matrix.
Expected utility for a=0 given z*: \(E[U(0,z)|z\le z*]\).
Expected utility for a=1 given z*: \(E[U(1,z)|z>z*]\).
Probability of acceptable offers: \(Prob(z>z*)\)
The user provides a method EUtility() which returns an array. The first element is a vector of expected utility differences. The second
is a vector of probabilities.
Approximate EV from a subsample of Θ using Keane-Wolpin (1994).
KW Approximation computes complete "brute force" \(\max\ v(\alpha)\) operator on
only a randomly chosen subsample of points in the endogenous state space \(\Theta\).
The results at these points are used to predict (extrapolate) choice probabilities
at the non-sampled states without \(\max v(\alpha)\) computations.
KW is useful when the endogenous state space Θ and the fully exogenous state space are large.
Recall that the exogenous vector ε contains discrete states that are IID and whose values
do not affect the transition of other variables directly (only indirectly through the choice α).
Semi-endogenous states, η, are not supported in this method.
A fatal error is produced if they appear in the model when a new KeaneWolpin object is created.
Different feasible sets A(θ) are allowed, but ...
The feasible set must be the same size at each state at a give clock setting t.
A warning message about this is issued if more than one feasible set exists.
When the dimension of the exogenous space is large KW can accomplish two things:
Drastically reduce the computational cost of value function iteration. This is because applying the max operator
to each value of ε at a non-sampled θ is replaced by applying it to only one value of ε
Drastically reduce the storage required by the model. This is because utilities and choice probabilities are
only stored for the single value of ε at non-sampled states
KW Approximation works well if the extrapolation method is good at approximating the value function at non-sampled states for
a relatively small number of sampled states.
At non-sampled endogenous states the one exogenous vector for which max operator is applied is the median/mean ε That is:
Assuming that each element of ε is mean zero and symmetric and the
discrete points also map into symmetric actual values, then the median discrete value corresponds
to the mean and median value of the ε vector
For example, if each element of ε takes on 5 values (0…4) then the 2 value is the median
and will be used for the extrapolation.
For this reason, it makes sense to have exogenous state variables take on an odd number of values when using
KW approximation in DDP.
The key elements of the approximation
The points in &Thetaf; to subsample at each clock setting, I::t.
The sampling is controlled by SubSampleStates(), which takes 3 optional arguments. See its
documentation for an explanation. At the subsampled points in θ all the exogenous states
are iterated over to compute the full value function. This value is stored as the
explained value for the approximation.
The explanatory values are also stored, in the default, the choice-specific values at the
MEDIAN point in the exogenous vector ε and the max of them.
The type of approximation used
Keane and Wolpin's preferred approach is to run a regression at the sampled states.
The regression is run to explain the value at sampled points then applied to non-sampled
states to predict the value.
The Specification of the approximation
KW's preferred specification is the default, but it can be replaced by the user
(no help yet available on this). The default is to run a linear regression in the V-v(α) vector and the
square root of the vector
Details
Brute force value iteration with exogenous and endogenous state variables can be written
Emax(θ) ≡ ∑ ε Ρε V(ε,θ)
which is the expected value of the maximum
V(ε,θ) = max α∈θ.A v(α;ε,θ)
When ε and θ both include many states, visiting every value of ε for every value of θ can be
expensive.
KW Approximation computes Emax for a randomly selected subset of states at a given t:
ΘKW(t) ⊂ Θ(t).
where Θ(t) is the subset of Θ with states at time t.
Then it interpolates Emax with a function (such as a linear regression) of the values of v(α;ε,θ) at a
single exogenous vector, ε, including the non-linear transformation:
maxE(θ) ≡ max α∈A(θ) v(α;ε,θ)
When ε = E(ε) then this is indeed the max at the expected exogenous shock.
Then, for points θ ∉ ΘKW(t) it extrapolates the value based on only
maxE(θ) and v(α;ε,θ).
.
KW (1994) report the extrapolation is sufficiently accurate in their model so that the simulation bias in parmeter estimates based on it is quite small.
The built-in components of KeaneWolpin use the KW (1994) "linear and square root" regression specification.
However, these components can be replaced by specifications or interpolating functions of the user's choice.
To do this, derive a new class from KeaneWolpin and substitute for the virtual components.
The amount of computations saved is easily approximated.
Let ε.D, Θ.D and ΘKW.D
denote the cardinality of the sets (following the notations used elsewhere).
Then the proportion of total max() operations avoided is
(ε.D-1)(Θ.D-ΘKW.D) / Θ.D
Initialization
Create a random subsample of states for each \(t\), denoted \(\Theta_t^S.\) This sampling can be controlled so no sampling takes place for some \(t\) and minimum and maximum counts are enforced.
Iteration
Follow these steps at each \(t.\)
For all \(\theta \in \Theta_t\):
Store the vector of action values at a single
\(\epsilon_0\) (e.g. the median or mean vector),
\(v_0(\alpha,\theta) = v(\alpha;\epsilon_0,\theta)\) and \(V_0(\theta)=\max v_0(\alpha,\theta).\)
(This is "maxE" in Keane and Wolpin 1994.)
if \(\theta \in \Theta_t^S\), compute \(V\) (or Emax), averaging over all values of the IID state vector \(\epsilon\): $$V(\theta) = {\sum}_{\epsilon} \ \left[ {\max}_{\ \alpha\in A(\theta)\ }\ U\left(\alpha;\epsilon,\theta\right)
+ \delta E_{\alpha,\theta\,}V\left(\theta^\prime\right)\right]\ f\left(\epsilon\right).\nonumber$$
Approximate \(V\) on the subsample as a function of values computed at \(\epsilon_0\). The default is to run the regression:
$$\hat{V}-V_0 = X\beta_t.\nonumber$$
The default specification of the row of state-specific values is
$$X =\left(\matrix{\left(V_0-v_0\right) & \sqrt{V_0-v_0}}\right).\nonumber$$
That is, the difference between Emax and maxE is a non-linear function of the differences in action values at the median shock.
For \(\theta \not\in \Theta^S_t,\) compute \(v_0(\theta),\) extrapolate the approximation:
$$V(\theta) = \max\{ V_0(\theta), V_0(\theta)+X\hat\beta_t\}.\nonumber$$
Iterate on Bellman's equation then switch to N-K iteration.
This method works for Ergodic environments or at a stationary period of a non-ergodic clock.
Implementation does not always work. Needs to be improved.
Initialization
Initialize as in Bellman iteration. Set a threshold \(\epsilon_{NK}\) for switching to N-K iteration.
Iteration
Begin with Bellman iteration, inserting a check for \(\Delta_t \lt \epsilon_{NK}.\)
When this occurs. set SetP=TRUE and SetPtrans=TRUE to compute the state-to-state transition on each iteration.
Replace step c in the Bellman Iteration that computes \(V_0 = Emax(V_1)\) with:
Solve for cutoffs as a non-linear system of equations in a one-dimensional choice problem.
This algorithm requires MyModel to be derived from OneDimensionalChoice.
The single Action Variable is denoted d:
α = (d)
The single continuous state variable is denoted z:
which are cutoff or reservation values for the values of z, the one-dimensional continuous state variable.
The optimal value of the choice d, denoted d*, is then
d* = j ⇔ z ∈ (z*j-1,z*j)
Note that z*-1 ≡ -∞ and z*d.N ≡ +∞
a* = a iff za-1 < ω ≤ za.
z-1 ≡ -∞
za.N ≡ +∞
The user writes routines that return ...
Initialization
Categorize each \(\th\) as an element of \(\Theta_Z\) if reservation values are needed there. The default is yes unless the user provides a Continuous() method. If so, create storage for \(z^\star\) at \(\th.\)
Iteration
Follow these steps at each \(t.\)
At each \(\th \in \Theta_Z\) initialize \(z\) solve for \(z^\star\) based on the user-provided \(Uz(A(\th);z).\) When completed store \(z^\star\) at \(\th.\) Use EUtility() that provides \(F(z^\star)\) and \(E[U|A(\th),z^\star]\) to compute \(V(\th)\) based on \eqref{EVz}.
For \(\th\) not in \(\Theta_Z\) compute \(V(\th) = U(\al;\th) + \delta EV(\thp)\) as in Bellman iteration but with a single action available.
Solve EV as as a non-linear system in a stationary EVExAnte environment.
In an ergodic system EV(θ) = EV'(θ) where EV(θ) is
the result of applying Bellman's equation to V'(θ). This can be written
as a system of non-linear equations to find the root:
EV(θ) - EV'(θ) = 0.
Rather than iterating on Bellman's equation from some initial EV'(θ), this approach uses
Newton-Raphson root solving to find the solution. This can be much faster, but less certain of success, than Bellman
iteration,
especially when δ is near 1.
Object ot iterate on Bellman's Equation, to solve EV(θ) for all fixed and random effects.
Solution is carried out by Solve()
Comments:
EV stores the result for each reachable endogenous state.
Results are integrated over random effects; results across fixed effects are overwritten.
pandv contains choice probabilities at \(\theta\) when complete.
This routine simplifies basic reservation value solving. Simply call it after calling CreateSpaces().
Its useful for debugging and demonstration purposes because the user's code does not need to create
the solution method object and call solve.
This would be inefficient to use in any context when a solution method is applied repeatedly.
This routine simplifies basic solving using Bellman value function iteration.
CAUTION All problems defined by fixed and random effects variables are solved. However, the solutions
cannot be used ex post for prediction or simulation because only the last group's problem remains in memory.
Instead, a solution method object must be nested to handle multiple problems. If there is are no group
variables then ex post prediction and simulation can follow use of VISolve()
Note: All parameters are optional, so VISolve() works.
It creates the ValueIteration object, calls the Solve() routine and then deletes the object.
It's useful for debugging and demonstration purposes because the user's code does not need to create
the solution method object and call solve.
This cannot be used if the solution is nested within some other iterative procedure because the solution method
object must be created and kept
Parameters:
ToScreen
TRUE [default] means output is displayed to output.
aM
address to return matrix of values and choice probabilities0, do not save [default]
MaxChoiceIndex
FALSE: print choice probability vector [default]
TRUE: only print index of choice with max probability. Useful when the full action matrix is very large.
TrimTerminals
FALSE [default] TRUE: states marked Type>=TERMINAL are deleted from the output
TrimZeroChoice
FALSE [default] TRUE: states with no choice are deleted
either AllFixed or a specific fixed group \(\gamma_f\) to solve
Rgroups
if AllFixed then must be AllRand. Otherwise, specific \(\gamma_r\)
Example:
vi = ValueIteration();
vi -> Solve();
AMGSolve
AMGSolve
Run
Apply the method (default is Bellman equation) at a point \(\theta\).
Compute ThetaUtility().
Compute the value of actions, \(v(A(\theta),\theta) by calling ActVal() or the
replacement for the actual method
Call thetaEMax() or replacment to store the value in the scratch space for \)V(\theta)$.
Call GSolve::PostEmax() or replacement to carry out post emax tasks, CCP smoothing, e.g.
Solve
Interate over the state space apply the solution method.
This is not called by the user's code. It is called for each point in the group space \(\Gamma.\)
Its job is to iterate over \(\Theta.\)
Set the setPstar for whether \(P^\star\) should be computed or not.
Compute utility
Iterate over \(\theta\) applying Bellman's equation or other solution method if replaced.
Call any functions added to the PostGSolveHookTimes
CCP
bandwidth
static decl bandwidth [public]
CCP
CCP :: CCP ( data , bandwidth )
Collect observed choice frequencies from a dataset.
Initialize EV() as a system of non-linear equations.
GSolve
AuxiliaryRun
virtual GSolve :: AuxiliaryRun ( instate )
AuxRun
decl AuxRun [public]
Flag set if extra loop after converging,
e.g. when computing semi-closed-form derivatives.
MaxTrips
decl MaxTrips [public]
Fixed limit on number of iterations.
Run
virtual GSolve :: Run ( )
Apply the method (default is Bellman equation) at a point \(\theta\).
Compute ThetaUtility().
Compute the value of actions, \(v(A(\theta),\theta) by calling ActVal() or the
replacement for the actual method
Call thetaEMax() or replacment to store the value in the scratch space for \)V(\theta)$.
Call GSolve::PostEmax() or replacement to carry out post emax tasks, CCP smoothing, e.g.
RunSafe
decl RunSafe [public]
TRUE (default): exit if NaNs encountered during iteration
FALSE: exit with IterationFailed
Solve
virtual GSolve :: Solve ( instate )
Interate over the state space apply the solution method.
This is not called by the user's code. It is called for each point in the group space \(\Gamma.\)
Its job is to iterate over \(\Theta.\)
Set the setPstar for whether \(P^\star\) should be computed or not.
Compute utility
Iterate over \(\theta\) applying Bellman's equation or other solution method if replaced.
Call any functions added to the PostGSolveHookTimes
Tolerance on value function convergence in stationary
environments. Default=DefTolerance.
HMGSolve
HMGSolve
HMGSolve :: HMGSolve ( caller )
Q
static decl Q [public]
F×1 array of CCPs.
Run
HMGSolve :: Run ( )
Apply the method (default is Bellman equation) at a point \(\theta\).
Compute ThetaUtility().
Compute the value of actions, \(v(A(\theta),\theta) by calling ActVal() or the
replacement for the actual method
Call thetaEMax() or replacment to store the value in the scratch space for \)V(\theta)$.
Call GSolve::PostEmax() or replacement to carry out post emax tasks, CCP smoothing, e.g.
Solve
HMGSolve :: Solve ( state )
Interate over the state space apply the solution method.
This is not called by the user's code. It is called for each point in the group space \(\Gamma.\)
Its job is to iterate over \(\Theta.\)
Set the setPstar for whether \(P^\star\) should be computed or not.
Compute utility
Iterate over \(\theta\) applying Bellman's equation or other solution method if replaced.
Call any functions added to the PostGSolveHookTimes
tmpP
static decl tmpP [public]
HotzMiller
AMiter
HotzMiller :: AMiter ( mle )
AMstep
static decl AMstep [public]
data
decl data [public]
EmpiricalCCP
HotzMiller :: EmpiricalCCP ( indata , bandwidth )
HotzMiller
HotzMiller :: HotzMiller ( data , mysolve )
Create a Hotz-Miller solution method.
Parameters:
(optional)
indata Panel object F×1 array of Q mappings 0 [default], no data sent
(V-vv)'
The default is to run the regression:
$$\hat{V}-V_0 = X\beta_t.\nonumber$$
The default specification of the row of state-specific values is
$$X =\l(\matrix{\l(V_0-v_0\r) & \sqrt{V_0-v_0}}\r).\nonumber$$
That is, the difference between Emax and maxE is a non-linear function of the differences in action values at the median shock.
Xmat
decl Xmat [public]
X matrix
Method
DefTolerance
static const decl DefTolerance [public]
Default convergence tolerance on Bellman Iteration for stationary
environments = DIFF_EPS = 10^{-8}.
a method that does the iteration over \(|Theta\) (user code does not set this)
Rgroups
decl Rgroups [public]
Either r or AllRan to solve for all random effects.
Solve
virtual Method :: Solve ( Fgroups , Rgroups )
Carry out a solution method.
Parameters:
Fgroups
either AllFixed or a specific fixed group \(\gamma_f\) to solve
Rgroups
if AllFixed then must be AllRand. Otherwise, specific \(\gamma_r\)
Example:
vi = ValueIteration();
vi -> Solve();
ToggleIterate
Method :: ToggleIterate ( ToggleOnlyTrans )
Toggle whether to Iterate on Bellman's Equation.
Parameters:
ToggleOnlyTrans
TRUE [default] also toggle OnlyTransitions so choice
probabilities do not apply in likelihood calculation.
This is used by DataObjectives that are maximized using "2 Stage" iteration.
ToggleRunSafe
Method :: ToggleRunSafe ( )
Toggle whether to check for NaNs in value iteration.
UseCurrent (-2) do not change
Volume level for output,
vtoler UseCurrent (-2) do not change
convergence criteria for stationary problems (fixed points)
MaxTrips UseCurrent (-2) do not change
integer: set the total number of iterations to perform (stationary problems)
0 : unlimited trips (stopping determined by tolerance and convergence)
NormType set the norm() itype argument for Ox's norm(v,itype) function
only applies to stationary phases of clocks.
UseCurrent(-2) do not change (default=2,Euclidean)
0: infinity norm, maximum absolute value,
1: sum of absolute values,
2: square root of sum of squares,
p: p-th root of sum of absolute elements raised to the power p.
-1: minimum absolute value,
Apply the method (default is Bellman equation) at a point \(\theta\).
Compute ThetaUtility().
Compute the value of actions, \(v(A(\theta),\theta) by calling ActVal() or the
replacement for the actual method
Call thetaEMax() or replacment to store the value in the scratch space for \)V(\theta)$.
Call GSolve::PostEmax() or replacement to carry out post emax tasks, CCP smoothing, e.g.
SaSGSolve
SaSGSolve :: SaSGSolve ( caller )
Solve
Interate over the state space apply the solution method.
This is not called by the user's code. It is called for each point in the group space \(\Gamma.\)
Its job is to iterate over \(\Theta.\)
Set the setPstar for whether \(P^\star\) should be computed or not.
Compute utility
Iterate over \(\theta\) applying Bellman's equation or other solution method if replaced.
Call any functions added to the PostGSolveHookTimes
Solve Bellman's Equation using brute force iteration over the state space.
Parameters:
Fgroups
DoAll, loop over fixed groups non-negative integer, solve only that fixed group index
Rgroups
MaxTrips
0, iterate until convergence positive integer, max number of iterations -1 (ResetValue), reset to 0.
Returns:
TRUE if all solutions succeed; FALSE if any fail.
This method carries out Bellman's iteration on the user-defined problem. It uses the ClockType of
the problem to determine whether it needs to find a fixed point or can simply work backwards in
time.
If UpdateTime[OnlyOnce] is TRUE (see UpdateTimes), then transitions and variables are updated here.
EV stores the result for each reachable endogenous state.
Results are integrated over random effects, but results across fixed effects are overwritten.
Choice probabilities are stored in pandv
ValueIteration
ValueIteration :: ValueIteration ( myGSolve )
Create a new "brute force" Bellman iteration method.
Parameters:
myGSolve
0 (default), built in task will be used. GSolve-derived object to use for iterating over endogenous states. User-code does not provide this.
Derived methods may send a replacement for GSolve. Ordinary users would only send an argument if they had
developed their own solution method.
Search | 
Up Level |
Project home
| 
Index
| 
Class hierarchy ]
Methods.ox
Global functions, enumerations
RVSolve
KWstages
AguirregabiriaMira : HotzMiller : Method : FETask : GroupTask : Task : DP
AguirregabiriaMira
bandwidth
