0. USER BELLMAN CLASS:    LS | Exteme Value  | Bellman
1. CLOCK:                 3. Normal Finite Horizon Aging
2. STATE VARIABLES
              |eps   |eta |theta -clock        |gamma
                e    s21      M      t     t''     r      f
       s.N     15      1     40     40      1      1      1
3. SIZE OF SPACES
                       Number of Points
    Exogenous(Epsilon)               15
    Endogenous(Theta)                40
    Times                            40
    EV()Iterating                    40
    ChoiceProb.track               1600
    Total Untrimmed               24000
5. TRIMMING AND SUBSAMPLING OF THE ENDOGENOUS STATE SPACE
                           N
    TotalReachable       820
         Terminal          0
     Approximated          0

C. Code Utility and Other Functions

To match the model specification in \eqref{LS1}, and in anticipation of empirical applications, earnings and utility are coded as separate functions:

    LS::Earn() {
       decl x;
1.     x = 1 ~ CV(M) ~ sqr(CV(M)) ~ AV(e);
2.     return exp( x*CV(beta) ) ;
       }
    LS::Utility() {
3.     return  CV(m)*(Earn()-pi) + pi;
       }

D-E. Solve and Use the Solution

The output of the prediction listed in A1 shows the agent works with probability 0.3982 in the first period. This integrates over the discrete distribution of earnings shocks and the continuous extreme-value choice smoothing shocks as well as optimal decisions. In the last period of life a large sample of (homogeneous) people would work 13% of the time and will have accumulated 7.52 years of experience on average.

The definition of the labor supply model corresponds roughly 1-to-1 with user code in niqlow. There has been no previous attempt to embed empirical discrete dynamic programming in a higher-level coding environment for even a simple class of models. It is true that one-time code for what has been shown so far is not complicated. Extensions of the model are shown which can be implemented with one or two lines in niqlow that would otherwise involve rewriting of the one-time code. Before discussing them, consider a side benefit of using a platform rather than purpose-built code: efficiency.

Efficient Computing

The labor supply model is a small problem, and a novice coder can implement it with straightforward code. However, efficient code is not necessarily simple or intuitive to write. As a novice builds on the small problem inefficiencies in their code can invoke the curse of dimensionality before necessary. When this happens they need to discover the inefficiencies and rewrite their code. These inflections points, where a novice's progress slows down while rewriting code, may determine when a project stops. That is, the point can be reached where the marginal cost of finding and eliminating inefficiencies exceeds the marginal value of additional complexity.

This section discusses three ways efficiency is automatically accounted for in niqlow. Many other strategies to increase efficiency, and to balance storage and processing requirements are encoded in niqlow. Models developed with niqlow do not avoid the curse of dimensionality, but they can delay it.

Time (and Memory) is of the Essence

A solution method could always assume that the DP problem is ergodic. The reason for not doing this is practical: past states would enter calculations unnecessarily. A practical DP platform exploits the reduced storage and computation implied by a non-stationary clock. It also must exploit methods to find fixed points quickly in stationary models.

As seen in the labor supply code (Line D.1), the clock is set by SetClock(). The clock controls how Bellman's iteration proceeds. If $t$ is a stationary phase (so it is possible that $t^\prime=t$) then a fixed point condition must be checked before allowing time to move back to $t-1.$ If, on the other hand, $t$ is a non-stationary phase then no fixed point criterion must be satisfied.

Further, empirical dynamic programming involves multiple stages which process (span) the state space, not just Bellman iteration. For example, once the value function has been computed the model can be used for simulation, prediction or estimation. These processes involve all values of time whereas Bellman's iteration only involves "today" and possible states "tomorrow". This creates another complication. While iterating on Bellman's equation the transition $P\st$ should map into only the possible next time periods. But when simulating outcomes the transitions must relate to model time.

niqlow addresses all these issues for the user without re-coding. It accounts for differences in how backward iteration proceeds and what future values are required to compute current values. It also uses different linear mappings from state values into points of the state space depending on whether it is accessing the value function or tracking model time.

Not all State Variables are Equally Endogenous

Recall that $P(\thp;\al,\th)$ is the primitive transition function. In niqlow the full transition is built up from the transitions of the state variables added to the model. If, after conditioning on $\al$ and $\th,$ a state variable $s$ evolves independently of all other variables its transition can be written $P(s';\al,\th).$ In niqlow $s$ is said to be autonomous. When all state variables are autonomous the transition is the product of the individual transitions: $$P^B(\thp;\al,\th) = {\prod}_{s\in\th} P(s';\al,\th).\label{PB}$$ In this baseline form each current state has a transition matrix associated with it: rows correspond to actions and columns to next states. The elements are the transition probabilities.

The baseline form may be either too simple or too complex than needed. On the one hand, state variables may not be autonomous. Two (or more) state variables are not autonomous if, after conditioning on the full state $\th$ and action $\al$, their innovations are still correlated with each other. In niqlow they must be placed inside a StateBlock which is then autonomous.

On the other, many state variables have simpler transitions that require less storage and recomputing than the baseline. The labor supply shock $e$ is IID. Its values next period is independent of everything including its own current value. There is no need to account for its distribution separately at each $(\al,\th)$ combination. If $e$ can be handled separately this removes 15 columns from the baseline transition matrix at each point $\th$.

The transition for experience $M$ does depend on both its current value and the action. Its transition cannot be factored out across states, but it depends on the current value of $e$ only through $\al$. If $e$ is isolated from $\th$ then fewer state-specific matrices are required to represent the transition $P^B$ in \eqref{PB}. This is why $e$ was added to the "exogenous" vector in D.7. This vector is denoted $\epsilon$ in the model summary to distinguish it from $\th$.

2 illustrates the full set of options for classifying state variables in order to economize on storage and calculations. It begins at the top where all state variables start out as generic except any co-evolving variables have been placed in a block represented by a single $s_i$. Each is filtered (by the user) into one of five vectors based on its transition. The two leftmost vectors, $\eps$ and $\eta$, contain variables whose transitions can be written $f(s^\prime_i)$ because they are IID. The two rightmost vectors, $\gamma_r$ and $\gamma_f$, contain variables that are fixed for the agent because they do not transit at all: $s^\prime_i=s_i.$ These vectors represent different DP problems not evolving states for a given problem. Examples are discussed later.

---

---

The middle category in 2 is the endogenous vector $\th.$ It contains state variables that are neither transitory nor permanent. Transitions for variables in $\th$ take the form $f(q ; \alpha, \eta, \th, \gamma_r, \gamma_f).$ That is, the transition can depend on the action and all other state variables except those in $\eps.$ If all state variables are placed in $\th$ then the result is the baseline transition \eqref{PB}. Naive code treats all state variables generically and computes the baseline transition at the cost of wasted computation or storage.

By definition, a variable placed in $\eps$ cannot directly affect other transitions. If an IID state variable directly affects the transitions of other state variables it can only be placed in $\eta.$ That is, $\eps$ and $\zeta$ vectors satisfy conditional independence but $\eta$ and $\th$ do not. An example is explained in Appendix B.

Extended Notation for DP Models

That is, $U()$ potentially depends on everything except $\zeta$ which only affects $v_\zeta().$ The primitive transition depends on everything except $\zeta$ and $\eps,$ so the same is true of the expectations operator over future states. The new categories of state variables allow niqlow to economize on storage and computing in large-scale DP projects. By the same token, these extra vectors can be ignored when solving a small DP.

At each point $\th$ the transition $P(\thp;\al,\th)$ must be stored. How big each of these matrices are depends on the transitions. niqlow updates these matrices before starting to solve the DP model so that the transitions can depend on parameters that are set outside the model.

The minimum additional required information that has to be stored at $\th$ is a matrix of dimension $\#A(\th) \times \#\eta$, where $\#\eta$ means the number of distinct values of the semi-exogenous vector. This space holds the value of actions $v\l(\al;\eta,\th\r).$ The exogenous vector $\eps$ is summed out at each value of $\eta$ which uses temporary storage for utility $U().$ Once $V(\th)$ is finalized, the matrix storing $v\l(\r)$ can be rewritten with the conditional choice probability matrix $P^\star\l(\al;\eta,\th\r)$ which has the same dimensions. This re-use of the same matrix cuts storage in half.

The group vector $\gamma$ accounts for multiple problems within the same model. This requires looping over the state space $\Theta$ for each separate problem. In naive code additional groups might expand the state space. However, niqlow re-uses $\Theta$ for each value of $\gamma$ which can drastically cut storage compared to code that stores all DP problems simultaneously. In niqlow almost no other information about the DP problem is duplicated at each point $\th$.

Not All States Are Reachable

Dedicated code for spanning the state space for the labor would account for unreachable states by changing the main loop, similar to this pseudo-code:

for ( t=39; t>=0; --t ) {
    for (M=0; M<=t; ++M) {
         ⁞
        }
    }

niqlow accounts for unreachable states for many state variable classes when CreateSpaces() is called to build $\Theta.$ For example, if the model has a finite horizon clock and includes an action counter like $M$, then only points that satisfy $M \le t$ are created. The user can override this if, for example, $M=0$ is not the initial condition. This one-time cost of deciding which states are reachable reduces storage requirements and lowers the ongoing cost of each state space iteration that may occur thousands of times during estimation.

Adding up Inefficiencies in Naive Code

First, a naive state space $\Theta$ would include $40\times 40 \times 15 = 24,000$ states. However, niqlow would average over the 15 values of the IID earnings shocks and store only a single value at each $\th.$ Next, niqlow reduces $\Theta$ to $40*41/2 = 820$ states through trimming of unreachable states. And it would only store the value function for $2\times 40=80$ points while iterating (one vector for $t+1$ and one for $t$). A naive solution might store $96,000$ numbers for action values and choice probabilities (whether reachable or not). Meanwhile, niqlow would store $820\times 15\times 2 = 24,600$ values, overwriting $\vv$ with $\Pstar.$

The baseline transition \eqref{PB} would (naively) require $820 \times 15= 12,300$ matrices, either to be stored or computed dynamically. Each matrix would be of dimension $2\times (15\times 40).$ With sorting into different state vectors only 820 matrices are required. The transition for $e$ is a single $1\times 15$ vector which is combined with the state-specific matrix when needed. Further, niqlow determines that only 2 states at most are feasible next period, so a $2\times 2$ matrix is stored.

Extensions

Adding or Modifying State Variables

First, suppose the earnings shocks should change from IID to correlated over time: $$e_{t+1} = 0.8e_t + z_{t+1}.\label{Tauch}$$ If the model were hard-coded in loops, this change would require a major rewrite. In niqlow it requires two simple changes to Build(). First, replace Zvariable() on line 4 with

4*.   e = new Tauchen("e",15,3.0,<0.0;1.0;0.8>);

Adding Variables That Create New Problems

In a second extension of the labor supply model the user accounts for differences across observed demographic groups and 5 levels of unobserved skill. One way to express this is to write the intercept term for earnings as a function of the agent's fixed characteristics: $$\beta_0 = \gamma_0 + \gamma_1 Hisp + \gamma_2 Black + \gamma_3 Fem+\gamma_4 k.\label{Beta0}$$ Now there are 30 dynamic programming problems, one for each combination of fixed factors that shift the intercept in earnings. Three more lines of code in the build segment expands the model for this change:

    LSext::Build() {
1.      Initialize(1.0,new LSext());
2.      LS::Build();
        ⁞
3.      x = new Regressors({"female","race"},<2,3>);
4.      k = new NormalRandomEffect ("skill", 5);
5.      GroupVariables(skill,x);
        ⁞
        CreateSpaces();
    }

This is why LS::Build() did not include calls to Initialize() and CreateSpaces(). They were placed in C. Now on line 2 it can be reused in D to set up the shared elements.

The Regressors class on line 3 holds a list of objects that act like a vector and can be used in regression-like equations such as earnings. The columns can be given labels and the number of distinct values are provided (in this case 2 and 3, respectively). The NormalRandomEffect on line 4 is like Zvariable() except it is a permanent value rather than an IID shock. To finish this extension another vector would be created for $\gamma$ and Earn() would be modified accordingly. The user might also modify utility to account for preference differences across groups as well.

Solution Methods

Bellman Iteration

Bellman iteration itself depends on details of the model, most notably the model's clock. 2 describes the algorithm and how allows properties of the clock to determine the calculations.

The form of the Emax operator itself also depends on the smoothing method, related to the presence of $\zeta$ in choice values. If no smoothing terms are present, Emax is simply the maximum of $v\l(\al;\eps,\th\r)$ over $\al\in A(\th).$ In general, the solution method relies on code related to the base class the DP model to handle it.

In purpose-built code, the simplest approach to working backwards in $t$ and spanning $\Theta_t$ involves nested loops over all state variables. In a language such as FORTRAN the depth of the nest would depend on the number of state variables in $ \th.$ Adding or dropping states requires inserting or deleting a loop.

A novice researcher may start with nested loops then realize there is an alternative. The discrete state space is converted to a large one-dimensional space with a mapping from the index back to the values of the state variables. This transformation is not trivial to modify as state variables are added or dropped from the model or when switching to algorithms that require different types of passes (such as the Keane-Wolpin approximation discussed below).

niqlow relies on a fixed depth of nesting independent of the length of vectors. The difference in the update stage is handled by a method of the clock. One segment of code works for all types of clock. Further, the same code handles tasks other than Bellman iteration. Each task, such as computing predictions, is a derived class with its own function that carries out the inner work at each state. New methods can be implemented without duplicating loops in different parts of the code. This approach has a fixed cost to create the state space each time the program begins. Hard coding loops over a pre-defined state space has not fixed cost at execution time but a larger cost of time spent re-coding as the model changes.

Variations on Value Iteration

Most methods attempt to reduce calculations relative to brute force methods. Rust (1987) used Newton-Kantorovich (NK) iteration, summarized in A1. This strategy applies to a model with an ergodic clock so the fixed point can be expressed as the root of a system of equations. It relies on the state-to-state transition defined in \eqref{Pqtoq}. First, initial Bellman iterations reduce $\Delta_t$, defined in 2, below a threshold. Then NK starts to update according to a Newton-Raphson step. Since NK is a class derived from ValueIteration it inherits the ordinary method but also contains the code to switch.

Hotz and Miller (1993) introduced the use of external data to guide the solution method as described in A2. This class is often referred to as CCP (conditional choice probability) methods, because it uses observed choices to obtain estimates of CCP and then map these back to values of $V(\theta)$ in one step without Bellman iteration. Aguirregabiria and Mira (2002) extends the one-step Hotz-Miller as summarized in A8. This effectively swaps the nesting introduced in Wolpin (1984): likelihood maximization changes utility parameters which are fed to $P^\star$ and then updates $V$ through Hotz-Miller.

The Keane and Wolpin (1994) approximation method, described in A3, is also derived from the basic brute force algorithm. It splits iteration over the state space into two stages. At the first stage it visits a subsample of states to compute Emax defined in \eqref{BE}. Information is collected to approximate the value function on the sample. At the second stage the remaining states are visited to extrapolate $V(\th)$ from the first stage approximation using information that is much less costly to compute than the full Emax operation.

A user coding the Keane-Wolpin algorithm from scratch faces extensive changes to all nested loops. As Keane and Wolpin (1994) report, the approximation can save substantial computational time but is not particularly accurate. It can be useful to get a first set of estimated parameters more quickly and then either increase the subampling proportion or simply revert to brute force.

A coder would likely copy their brute force loops and then "hack" them to carry out the two Keane-Wolpin stages. Now the code has two nested loops that need to be kept in synch as the model changes. In niqlow the two algorithms can be compared by simply adding three lines of code regardless of other changes in the model:

    ⁞
     vi = new ValueIteration();
     kw = new KeaneWolpin();
1    vi->Solve();
2    SubSampleStates ( 0.1, 30 , 200 );
3    kw -> Solve();

Another common problem combines Bellman iteration with the calculation of reservation values of a continuous variable $Z,$ where $Z \sim G(z)$ and is IID over time. Like the choice-specific value shocks $\zeta$, values of $Z$ are neither stored nor represented as an element of the ordinary state vector. Instead, $G(z)$ affects equations in the solution method and ultimately choice probabilities. The basic reservation value method is described in A4 along with the few changes required to convert the labor supply model to a reservation value problem.

Parameters, Data, and Estimation

Overview

Most empirical DP uses variants of the nested algorithm introduced in Wolpin (1984) and illustrated as a two-sided feedback loop in 3. Given $\psi$ the DP model is solved and outcomes (CCPs) produced. Parameters are changed by a numerical optimization routine. It is this feedback that MaCurdy (1981) avoided by approximating the Lagrangian rather than computing its value based on the current value of the regression coefficients.

4 illustrates more layers of dependency. The top (or outer) level is an optimization algorithm that controls $\psi.$ The parameters are tied to an econometric objective at the next level down. Levels 1 and 2 correspond to the two sides of 3 which hides the layers below.

The objective relies on a data set (level 3) which must be organized to match the model to external data allowing for issues such as unobserved states and measurement error. Model outcomes and predictions use the DP solution method (level 4). Finally, the base of the pyramid is one or more DP models.

Each level of 4 must interface with the adjacent levels. The DP model must also access the structural vector $\psi$ controlled at the top level. niqlow provides classes for each level and the connective tissue between them. This integrated approach then supports automated construction of estimation problems. That is, the user need not write code for higher levels and can focus on the code at the base.

---

---

Estimation methods such as Aguirregabiria and Mira (2002)'s pseudo MLE and Imai et al (2009)'s MCMC estimation avoid computations within stages and swap the order of levels in 4. Dedicated purpose-built code to implement them looks very different than nested solution code. However, in niqlow the levels of 4 are modular. They can be re-ordered as long as different connections between objects representing each level have been written.

Previous sections used the labor supply model to demonstrate how to build up a DP as objects derived from the Bellman class, which corresponds to Level 5 of the pyramid. DP solution methods represent level 4. This section discusses how Levels 2 and 3 are represented in niqlow.

---

---

Outcomes and Predictions

Empirical work requires outcomes to be defined that mask components of the full outcome. Outcomes must also be put in a sequence to create the path of an individual agent. Unlike reduced-form econometrics empirical DP data cannot generically be represented as a matrix or even a multi-dimensional array of numbers. Instead, in niqlow they are stored as linked-lists of objects of the Outcome class. The next section builds on outcomes to define "generic" likelihood functions. The niqlow code is shown to estimate the simple labor supply model without the user coding the likelihood.

A prediction, on the other hand, is the expected value of outcomes conditional on some information. The basic prediction integrates over all contemporaneous randomness conditional on $\th$ and the permanent variables: $$E\l[Y\l(\th;\gamma_r,\gamma_f\r)\r] \equiv {\sum}_\eta {\sum}_\eps {\sum}_{\al\in A(\th)} f_\eta(\eta) f_\eps(\eps) P^{\star}(\al;\eps,\eta,\th,\gamma_r,\gamma_f) Y(\psi).\label{EYF}$$ This is what an agent expects to happen (in a vector sense) given that $\th$ has been realized but the IID state variables have not. The continuous shock $\zeta$ is integrated out by the choice probability $P^{\star}.$

As with outcomes, the information not available in the external data will not always match up with the conditional information in \eqref{EYF}. And predictions must be put in sequence to create a path. After discussing likelihod 5.5 returns to prediction and how they are used to construct GMM objectives.

Likelihood

niqlow can determine which likelihood to apply automatically from the data read into an OutcomeDataSet object. The appropriate formula is computed without further user coding. Since each level of information relaxes the previous one, the PO algorithm could be applied to the other forms. Using the more restricted forms when possible speeds computation.

F: Full Likelihood

To build a likelihood function using full information, first assume there are no random effects variables so it is irrelevant that $\gamma_r$ is missing in $Y^F$. Then the likelihood for a single single observed outcome is the probability that the observed action would be taken: $$L^F\left(\ Y(\hat\psi)\ ;\ \hat{Y}\right) = P^{\,\star}(\hat\al; \hat\eps,\hat\eta,\hat\th,\cdot,\hat\gamma_f ).\label{LF}$$

This compact expression requires some explanation. The superscript has shifted to $L^F$ and the model outcome is shown with no superscript to avoid duplicate notation. On the right hand side is the conditional choice probability. This takes one of the forms in \eqref{CCP1}-\eqref{CCP3}.

Since the observed outcome includes the action and all state vectors, their data values are inserted into the conditioning values. niqlow integrates data and predictions so it inserts the observed values for the user. The model and external versions of the data illustrated in Level 3 of 4 correspond to $Y(\psi)$ and $\hat{Y},$ respectively. Calculations such as \eqref{LF} then enter an overall objective at Level 4.

All But IID Likelihood

This leads to the third outcome type: $$Y^{IID}(\psi) \equiv \l(\matrix{y& \al & \cdot &\cdot &\eta&\th &\cdot &\gamma_f}\r).\label{YIID}$$ The auxiliary outcomes have been added to the front of the outcome. The agent has more information than $y$ so it is unnecessary for them to condition choices or transitions on $y.$ On the other hand, when the data is less than the agent's defined in \eqref{Yagent}, $y$ can contains additional information.

Both the smoothing shock $\zeta$ and the IID state $\eps$ affect the DP transition only through the choice of $\al.$ This means the IID outcome is sufficient to predict the next outcome as the agent would, namely using $P\st.$ Likelihood of a sequence of individual outcomes can be constructed that integrates only over the IID elements. The other IID vector in \eqref{ExtNotation}, $\eta,$ does not satisfy this condition because it can directly affect the transitions of other state variables.

The likelihood of an outcome when $\eps$ is unobserved is an expectation of $L^F:$ $$L^{IID}(\ Y(\hat\psi)\ ,\ \hat{Y}\ ) = {\sum}_{\eps} f\l( \eps\, |\, y=\hat{y} \r) L^F\l(Y,\hat{Y}\r). \label{LIID1}$$

The choice probability component is weighted by the conditional probability of $\eps$ given that it generates the observed value of $y.$ In this case \eqref{LIID1} can be discontinuous in $\psi.$ In the case of the labor supply model only observed earnings on one of the 15 points of support for $e$ would have positive likelihood. Further, if the set of discrete values of $\eps$ consistent with the model changes with a small change in $\psi$ it causes a jump in likelihood.

It is standard to add measurement error to observed auxiliary values (and other states or actions) in order to smooth out the IID likelihood. The measurement error is ex post to the agent's problem so it enters only at this stage: $$L\left(\ Y(\hat\psi) \ ;\ \hat{Y}\ \right)=P^{\,\star}(\hat{\al} ; \cdots )f(\eps)L_y(\hat{y},y). \label{LIID}$$

Here $L_y()$ is the likelihood contribution for the data given the model's predicted value for the auxiliary vector $y.$ As shown below, niqlow includes built-in classes to add normal linear or log-linear noise to any outcome's likelihood contribution derived from the Noisy class.

Path and Panel Likelihood

Let $g(\gamma_r)$ be the probability distribution of the random effects vector. The distribution can depend implicitly on $\gamma_f$ and $\psi.$ Multiple paths with the same fixed effects are stored in a "fixed panel." These fixed panels are concatenated across $\gamma_f$ in a "panel." A panel might hold simulated data only, but if external data will be read in then it is represented in niqlow by the OutcomeDataSet class.

Suppose the information available at a single decision is either full (F) or everything-but-IID as defined above. Let $\tau \in \{F,IID\}$ indicate which type of data are observed. Then the likelihood for a single agent's path is $$\eqalign{ L\l(\ \{Y\} \l(\hat\psi\r) \, , \,\{\hat{Y}\}\ \r) &= \cr {\sum}_{\gamma_r} g\l(\gamma_r \r){\prod}_{s=0}^{\hat{T}}& \l\{ L^{\tau}(Y_s,\hat{Y}_s)\times \left[ P\left(\hat{\th}_{s+1}; \hat{\al}_s,\hat{\eta}_s,\hat{\th}_s\right)\right]^{I\{s\lt \hat{T}\}}\r\}. \cr}\label{Lpath}$$

The rightmost term is the model probability for observed state-to-state transitions which only applies before the last observation on the path.

Efficient computation and storage of both the DP solution and this likelihood requires coordination. Placement of the "nested" solution algorithm must be exact. When there are fixed effects the model must be solved for the $G$ combinations of $\gamma_r$ and $\gamma_f.$ As discussed earlier, the endogenous state space $\Theta$ is not duplicated for each combination of permanent values. Otherwise storage requirements would multiply $G$-fold. On other hand, this means each combination must be fully processed for a given structural vector $\psi$ before proceeding to the next group. The fixed panel must initialize a vector to contain $L()$ for each of its members and then store partial calculations for $\gamma_r.$ Only when the outer summation in \eqref{Lpath} is complete can the log of the path likelihood be taken and summed across paths to form the log-likelihood. This must then be repeated for each fixed effects vector.

2-Stage Estimation

More generally, in two-stage estimation each parameter is given one of the three markers by the user. The full parameter vector then contains three lists (sub-vectors): $$\psi = \left(\matrix{\psi_0& \psi_p & \psi_u }\right).\label{Subvecs}$$

The first, $\psi_0,$ contains parameters held fixed throughout estimation, including weakly identified parameters or ones set through "calibration." The second vector, $\psi_p,$ includes parameters that affect only transitions. The final vector, $\psi_u,$ contains parameters that affect utility and the discount factor $\delta$ if it is estimated.

At the first stage only $\psi_p$ is estimated using the limited information path likelihood: $$L^{1st} (\ Y(\hat\psi)\,,\,\hat{Y}\ ) = {\prod}_{s=0}^{\hat{T}} 1\times P(\hat{\th}_{s+1}; \hat{\al}_s,\hat{\th}_s)^{I\{s\lt \hat{T}\}}. \label{FirstStage}$$

This is the same as \eqref{Lpath} except "1" replaces the model's generated CCP. The likelihood conditions on the observed choice can be computed without solving Bellman's equation for $P^{\star}.$ The rest of the model setup is still required to compute the outcome-to-next-state transitions along the path. After estimating $\psi_p$ they are fixed at those values and $\psi_u$ is made variable. The likelihood reverts to \eqref{Lpath}. In niqlow a user implements 2-stage estimation with a few lines of code.

Type PO: Partial Observability

More generally, partial observability of the outcome creates a difficulty for computing the path likelihood forward in time. Let $q$ be a state variable in the endogenous vector $\th$ that is unobserved in the data at outcome $Y_s.$ It could be missing systematically as a hidden state or incidentally for this outcome alone. It has a distribution conditional on past outcomes that can be computed moving forward. So the contribution to likelihood at $s$ can be computed by summing over possible values of $q$ weighted by its conditional distribution. However, the distribution at $s+1$ depends on the distribution of $q$ at $s.$ So the distribution must be carried forward in the calculation. Each unobserved value creates more discrete distributions to sum over.

Even if all the unobserved states are tracked and their joint distribution across missing information computed, this can still be inadequate to compute the correct likelihood. With $q$ missing at $Y_s,$ suppose also that the future point $s+k$ it equals $q^{\star\star}$ in the data. However, that value happens to have zero probability of occurring if $q=q^\star$ at $s.$ For example, if a stock was only at $q^\star$ at $s$ it may be impossible for the stock to grow to $q^{\star\star}$ by $s+k.$

Since $q$ was unobserved at $s$ the forward likelihood included $q^\star$ in the sum and the distribution. Now the contributions of paths that start at $q^\star$ must be eliminated between $s$ and $s+k.$ The likelihood must sum over only unobserved states conditional on past observed outcomes and future consistent outcomes. In general there is no way to calculate the likelihood forward when endogenous states or actions are unobserved.

Ferrall (2003) proposed 6 to handle partial observability. Instead of computing the likelihood of a path going forward in time, it is computed backwards starting at the last observed outcome ${\hat T}.$ At any point $s$ the likelihood computed so far is not a single number. Instead, likelihood is a number attached to every state at $s$ that is consistent with the observed data from $s$ forward. This avoids the trap of following paths forward that end up inconsistent with later outcomes.

A one-dimensional vector of likelihoods, denoted $L^{PO}_s,$ contains the information to track likelihood contributions farther in the future. If, at a particular $s$, the IID-level outcome is observed then this vector collapses to a single point. Encountering missing values for a smaller $s$ will expand the $L^{PO}_s$ to a vector again. The algorithm works for the special cases of Full or IID observability as well, but moving backwards is slower than algorithms that can move forward in $s$. When data are read into niqlow the user flags which variables are observed, and all other variables are treated as unobserved. As the data on the observed variables is read in incidental missing values are also detected. Then each observed path can be assigned one of the three tags $\tau \in \{F,IID,PO\}$ and the most efficient path likelihood is then computed.

Estimating the Labor Supply Model

All work decisions are observed and there are no gaps in model time ($t_{s+1}=t_s$). Actual earnings, $E^a,$ depend on the observed work choice: $$E^a = \cases{ E() & if m=1\cr . & if m=0.\cr}\nonumber$$ The earnings shock $e$ is unobserved when not working. It could be inferred from $E^a$ except observed earnings, $E^o,$ include log-linear measurement error: $$E^o = e^{\nu}E^a,\quad \nu\sim N(0,\sigma^2).\nonumber$$

The auxiliary vector $y$ appearing in \eqref{YIID} contains $E^o.$ The observed path does not necessarily start at $t=0$. If it does, then observed experience $M^o$ equals actual experience because it can be computed from past work decisions before sending the data to niqlow. Otherwise, we'll assume that $M^o_0$ equals initial experience (acquired perhaps through retrospect questions). In this case \eqref{Lpath} is the likelihood for one observation with type $\tau=IID.$

To bring the simple labor supply model to this data, derive from LS a new class:

class LSemp : LS {
 static decl obsearn, dta, lnlk, mle, vi;
 static ActualEarn();
 static Build();
 static Estimate();
    }

    LSemp::Build() {
1.      Initialize(1.0,new LSemp());
2.      LS::Build();
3.      obsearn = new Noisy(ActualEarn);
4.      AuxiliaryOutcomes(obsearn);
5.      CreateSpaces();
        }

The required Initialize() function is called on line 1. The only difference with the earlier call is that the state space $\Theta$ consists of LSemp objects not LS objects. The empirical version of the model shares the same set-up as the earlier version, so LS::Build() can be called on line 2. Observed earnings is created as an object on line 3. The Noisy class adds measurement error to the argument, actual earnings $E^a$ coded as a static function:

LSemp::ActualEarn() {
    return m->myEV() ? Earn() :.NaN;
    }

The function for actual earnings requires some explanation. It uses Ox's .NaN code for missing values when not working ($m=0$). This is determined by m->myEV() instead of, say, $CV(m),$ which was explained earlier as returning the current value of m. That is used during Bellman iteration when $m$ is taking on hypothetical values.

At the estimation stage, however, the DP model must use the external value of actions to provide a predicted value of earnings. This was seen in the likelihood \eqref{LF} where the observed action $\hat\al$ enters the choice probability. niqlow handles this through its Data-derived classes that set values shared with the DP model. myEV() is another method of action and state variables. It retrieves the realized value of the variable if called during a post-solution operation such as likelihood calculation.

In ordinary econometrics, such as panel IV techniques, an endogenous variable is its own value. Only the value in the data is relevant. The need for CV() and myEV() reflects the complexity of nested estimation algorithms handled by niqlow and illustrated in 4. At the solution stage variables must take on hypothetical values to solve the model. At the estimation stage values from the external data set must be inserted into the contingent solution values such as the CCP.

Estimate() contains the code to compute the MLE estimates of $\beta,$ and selected output from it is discussed in Appendix D.

    LSemp::Estimate() {
1.      beta = new Coefficients("B",beta);

2.      vi = new ValueIteration();

3.      dta = new OutcomeDataSet("data",vi);
4.      dta -> ObservedWithLabel(m,M,obsearn);
5.      dta -> Read("LS.dta");

6.      lnlk = new DataObjective("lnlk",dta,beta);

7.      mle  = new BHHH(lnlk);
8.      mle -> Iterate();
        }

In the basic model $\beta$ was an ordinary vector. Its value was set inside LS::Build() in D. In this version the goal is to estimate $\beta$ from data. So on line 1 its value is replaced. It now holds an object of the Coefficients class. This is designed to hold a vector of freely varying parameters like regression coefficients. A constant vector is needed as default for starting values. Since beta contains a vector before this line, its current value can be sent as the initial vector. By the end of line 1 that vector is stored internally in the object and beta now contains an object not a vector.

The code for the true earnings function LS::Earn() in E already included CV(beta). So when this model is estimated it will retrieve the values from the parameter object under the control of an optimization algorithm.

In the first use of the labor supply model the simple VISolve() function carried out Bellman iteration once. Its use of a solution method object was hidden from the user. Now that the solution method is nested within a likelihood calculation it is not enough. The vi member holds the value iteration object (line 2). Somehow this object must be embedded (nested) within the estimation procedure.

Lines 3-5 in K create the data set object. As illustrated in 4 this handles both external data and model predictions. The dataset object also needs to know which model outcomes are in the external data (implying other variables should be treated as unobserved). Line 4 says that $m$, $M$, and observed earnings are in the data and will have the same labels as their objects. Now the data can be read from a Stata file on line 5. A column in the external data must hold the ID of each path and the value of $t$ so that a panel of paths can be created. Since they were not set explicitly in the code, the default labels will be used and must appear in the data file to avoid an error. As Read() reads in the data it determines whether there are missing values along the path and what type of variable is missing (an element of $\eps$ versus other vectors that imply partial observability.) By line 6 the class of likelihood function for each path in the data has been determined.

The DP solution method created on line 2 of K has already been embedded at line 3, because vi was sent as the second argument when dta was created. Whenever a new likelihood evaluation is needed vi->Solve() will be called to re-solve the model. This is the "nesting" of the solution algorithm inside the sample likelihood.

Line 6 creates the econometric objective. DataObjective expects a data set such as dta to be sent to it. It has a built in method to compute the log likelihood. The objective is the "home" of the parameters to be estimated, $\beta.$

Lines 7 creates the algorithm to carry out the outer iteration of maximizing the likelihood. Note that the DP model is analogous to the econometric objective lnlk, and the Bellman solution method vi is analogous to the optimization algorithm. One difference is that the user's code must send lnlk to the optimization algorithm, whereas the value function method does not take an object.

This example uses the BHHH algorithm (line 7), which is Newton's method except the outer product of the likelihood's Jacobian matrix is used to approximate the Hessian. This happens without any special coding from the user because any Objective object has two methods for evaluating itself. One returns a scalar. The other returns a vector to be aggregated into the scalar. In this case, lnlk returns the vector of log-likelihoods for each observation. In turn BHHH uses the vector version to compute the matrix of partial derivatives with respect to $\psi$ and then the outer product.

Line 8 in K calls the Iterate() method of the BHHH object to maximize the log-likelihood starting at the initial values of $\beta.$ Values read in from a file can supplant the hard-coded starting values of beta. This is analogous to Method classes for different DP solution methods, each having a Solve() function to carry out the task on demand.

An additional 14 lines of executable code were used to estimate the simple labor supply model. The chains of interactions in the code segment above matches the stacking illustrated in 4. The DP model itself (the base level) is implicit, because only one state space can be defined. The user's code could substitute different classes at each level of the process without changing any other code. It can create different solution methods, load different data sets, or compare different optimization algorithms on the same objective. None of these changes require changing the underlying code for the DP problem as long as CV() and other preparations discussed above have been used to make the code ready for alterations.

Moments and Predictions

A prediction, denoted $E[Y],$ was defined in \eqref{EYF} as the integral over current IID randomness and conditional choices. Predictions always integrate over permanent random effects ($\gamma_r$) since they are treated as unobserved. Predictions always condition on permanent fixed effects ($\gamma_f)$ since they are treated as observed. Thus, for each $\gamma_f$, there is a single predicted path, although multiple outcome paths (used for MLE and data simulation) can be generated and/or observed in data for each $\gamma_f$. Concatenating prediction paths across different fixed effects produces a prediction panel.

Sequence of Predictions

Now consider the distribution of $\th$ after $s$ decisions, $Q_s(\th).$ This distribution starts from $Q_0()$ and then accounts for the distribution of actions at 0 and the subsequent states at $1$ and so forth until $s$ actions have been taken. The prediction is the expected outcome over optimal conditional choice probabilities: $$E_s[Y] = {\sum}_{\th\in \Theta}{\sum}_{\al\in A(\th)} P^{\star}(\al;\th)Q_s(\th) EY(\al,\th).\label{EY}$$ Computationally this involves a loop over the state space and an inner product of the outcome and conditional choice probability vectors. Within the same loop, the distribution over states in the next period can be computed recursively.

$Q_{s+1}$ is accumulated over current state transitions: $$Q_{s+1}\l(\thp\r) = Q_{s+1}\l(\thp\r) + {\sum}_{\al\in A(\th)} P\st Q_s\l(\th\r).\label{Qs}$$ Once the prediction and next stage's distribution are computed, $s$ is incremented and the process repeated to form a path of predicted outcomes of a desired length.

Unlike path likelihood, which must deal with missing information along the path, $Y$ can be treated as the full outcome. A1 displays one element of $E_s[Y]$ vector, namely the action $m,$ along the path prediction for the simple labor supply. The full vector would also include the state $M.$ The exogenous states would by definition have time-invariant predictions and are not recorded in $E[Y].$ However, they can enter auxiliary outcomes which would be included in the prediction.

The external data may have missing information. Continue to define ${\hat Y}_s$ as the data but now is not an individual path but an average across ex-ante identical agents sharing a value of $\gamma_f.$ It includes only observed components of the outcome. To match empirical moments (averages) to the model prediction $E[Y]$, a masking matrix $M_s$ selects columns of the full outcome that are observed at stage $s.$ Then the differences between the predicted moment and the empirical moment is $$\Delta Y_s = \l( {\hat Y}_s - M_s{\hat E}_s\l[Y\r]\r).\label{DeltaYs}$$

Masking is applied to the predicted moment only to make it conform to the external data vector.

GMM Objectives

To form an econometric objective the vector of differences is aggregated into a scalar value. There are 3 types of GMM objectives built into niqlow. The simplest case is a weighted sum of all differences: $$J_0 = -{\sum}_{s=0}^{\hat T} \l\| {\sum}_{j=1}^J \omega_{sj} \Delta Y_s^j \r\|.\label{J0}$$

The negative sign appears because objectives are maximized. The user specifies the weights $\omega_{sj}$ to place on each observed moment $j$ at each point $s,$ including an option to set all weights equal. These weights can also be ad hoc "importance" weights. Further, the number of observations averaged at each $s$ can be read in so that weights also account for precision in the empirical moments. If the model includes different observed groups ($\gamma_f$) there is a summation over their separate values of in \eqref{J0} and the subsequent forms.

The next option is to specify a matrix or matrices to account for contemporaneous correlations between observed moments: $$J_1 = -{\sum}_{s=0}^{\hat T} \| \Delta Y_s\ \Omega_s\ \Delta Y_s\|.\label{J1}$$

Finally, efficient GMM requires that a weighting matrix for the full path of moment differences: $$J_2 = \Delta\{Y\}\ \Omega\ \Delta\{Y\}.\label{J2}$$

The efficient matrix $\Omega$ can be computed from individual data if available and the empirical moments are computed from them. However, this is not helpful if GMM is used because only moments over individual paths are available. niqlow includes procedures that will simulate individual paths and then compute $\Omega$ from them using a first-stage (consistent but inefficient) parameter vector.

The previous section discussed all the code involved in estimating the labor supply model using MLE. Much of that code would remain if GMM were used instead. The key differences would be a single statement:

    momdta = new PredictionDataSet(UseLabel,"avgdata",vi);

Roadblocks

In Wolpin (1984) the binary choice prevents dynamics from concentrating into a Lagrange multiplier. Estimated parameters are not just coefficients on observables. Instead they enter choice probabilities through their effect on the value function. Choice probabilities are contingent, and the observed choices determine which ones enter the likelihood. When faced with these data and parameter interactions no software emerged to assemble a DP model and derived its empirical content from pre-defined components.

Empirical DP remains compute bound: processor speed is a limiting factor to the scale of the model. Custom-built compiled code can be the most efficient in execution speed. But the interlocking roles of data and parameters illustrated in 4 results in "hard-wired" details for the given model. A change in the model requires rewriting deep code.

Reduced-form estimation is not as computationally intensive and lacks hidden layers such as solving Bellman's equation to compute choice probabilities. This research could move to platforms such as Stata even when it was restricted in storage and speed compared to purpose-built compiled code. Meanwhile, empirical DP continued to require fast single-purpose code for each application. Once complete there is little incentive to return to the code and make it general, which would be difficult for the reasons given above. Even when code is freely shared with others it is of limited use to someone building a different but related model. So, unlike panel IV estimation represented by MaCurdy (1981), no natural shift occurred from bespoke PP code to use of off-the-shelf OOP components for the empirical DP approach introduced by Wolpin (1984).

To overcome this roadblock, niqlow did not start from purpose-built code. Instead, it started with OOP representations of generic DP models before solving a single specific model. Classes encode features of elements independent of the rest of the model so the model can be built from components. This creates computational overhead compared to, say, nested loops in a compiled language. However, as discussed in 3.3, efficiencies in the design of niqlow are coded once-and-for-all. Integrated data and optimization tools further reduce specialized coding.

Custom-built, expertly-tuned code for a single problem will always run quicker than equivalent code in niqlow if the starting line for the comparison is set to when programs start executing. If, instead, the starting line is set at the first attempt to build the model then time-to-completion is likely to even out.

Conclusion

There is little evidence from the literature that sharing purpose-built code for a published paper has promoted verification or direct extension of already published empirical DP work. Referees of empirical DP models are also rarely in a position to replicate results reported in manuscripts. Finding bugs in purpose-built code for a large empirical DP model is never going to be practical. At best, referees can point to discrepancies in output that might be caused by mistakes that the authors can work to resolve.

Results generated from empirical DP papers have played at best an indirect role in shaping policy. One reason for this is skepticism among researchers not performing empirical DP themselves. Results are opaque, not independently verified, and rarely subject to robustness checks. Thus, a conundrum exists: empirical DP results are not trusted in part because they are based on purpose-built code that is rarely if ever verified. Because the results are not trusted they rarely play a role in policy debates. Since they do not influence policy, unverified results stay unverified for lack of relevance.

Perhaps niqlow is a step towards closing the structural coding gap, which lowers the cost of producing new results and replicating old ones. This in turn may make results more relevant and comparable to reduced-form methods already produced by portable and easy to use code.

Notes

Appendix

Appendix A. Example of PP vs. OOP Coding

Most readers have probably coded an objective and then called a built-in optimization procedure to optimize that function. Marshall is a specialized version of that general problem.

The PP package documentation explains how users should code $U()$ in order to interact with tools in the package. The user codes u(x), and sends it to a built-in procedure of the form demand(u,p,m). That procedure uses algorithms to compute $x^\star(p,m).$ Suppose the user wants to use the Cobb-Douglas function, $U(x) = {\sum} \ln x_i.$ Using "pseudo-code" the key parts of the user's program might look like:

#uses Marshall
u(x) {
    return sum_i ln(x[i])
    }
qdemand = demand(u,prices,income)
print("x* = ",qdemand)

class Consumer {
    members
                xstar, p, m
    methods
                demand()
                budget(p,m)
        virtual u(x)
    }

The package comes with the Cobb-Douglas function set as the a default to demonstrate the package without any coding. Now $u()$ is coded as a method of the Consumer class:

Consumer::u(x) {
    return sum( log(x) )
    }

#uses Marshall
    ⁞
agent = new Consumer()
agent -> budget(prices,income)
agent -> demand()

The code so far uses a built-in utility. To use, say, a CES utility the user creates a class derived from Consumer.

class CES : Consumer {
    members  a
    methods  u(x) CES(ina)
    }
CES::CES(ina) {     a = ina    }
CES::u(x)     {     return sum_i (x[i]^a)^(1/a)    }
⁞
agent = new CES(0.2)
agent -> demand()
⁞

When demand() is invoked, the user-provided CES::u() is called instead of the default version written by the programmer. The user has not changed, and perhaps cannot even see, the original code for demand(). This is because the programmer marked Consumer::u() as virtual. By doing so the programmer gives the user a controlled ability to change the underlying code through insertion a replacement function. In the PP packaged this injection of code was accomplished by passing u to the demand function. When many functions need to be replaced in many algorithms the PP framework can become unwieldily and unreliable. The OOP approach scales more efficiently for both the programmer and the user with problem complexity. With OOP it is easier to ensure the right data and the right functions are being used within the package.

The programmer can create a taxonomy of classes for the user to choose from. In this simple case, Marshall might not just have a single Consumer class. It might have child classes for different classes of utility. The user can then start with one of those classes to specialize or extend it for their model. An OOP package can provide the user with a menu of options, an important part of the niqlow approach to empirical DP.

Appendix B: Code Segments and Additional Explanation of niqlow Features.

$CV()$ and $AV()$ in Code Segment E

Recall that utility is treated as a vector valued function corresponding to the feasible set $A(\th).$ Since $m$ is an action variable its current value is not a scalar at $\theta.$ In this simple one-choice model the current value is always the same: $CV(m) = \l(\matrix{0\cr 1}\r).$ If other actions were added, or if constraints on feasible actions were imposed on choice, the current value of $m$ would be a different length because the number of distinct actions $\alpha$ would change.

Since $M$ is a simple count variable, it takes on values like a loop counter or index. The earnings shock $e$, is also like a loop counter, so its Its counter value ranges from 0 to 14. However, $e$ is a discretized normal random variable and the counting values are associated with both positive and negative real numbers, i.e. quantiles of the standard normal distribution, such as $-1.282,$ the 10th percentile of $N(0,1)$. The user's code can carry out these transformations of the integer value e.v, but niqlow can track actual values of variables for the user.

The actual values of an object are stored as vector member actual. Thus the actual value at any point is actual[v]. The current value is an index into the actual vector. The function AV() function retrieves this value, so when $CV(e)=3$ $AV(e)$ might equal -1.282. For $M$ the actual vector is (0 1 … 39) and actual[v]=v. That is, the default is that AV(s)=CV(s). Only in the case like a discretized normal will there be a difference. The user can set actual values for their state and action values and can make them dependent on structural parameters.

Semi-Exogenous State Variables

acc = new BinaryChoice("acc");
offer = new  LogNormal("offer" , 10);
curw = RetainMatch(offer,acc,1,0);
ActionVariable(acc);
SemiExogenousStates(offer);
EndogenousStates(curw);

Additional Actions and Restricted Choices

LSext::Build() {
    ⁞
    LS::Build();
    s = new BinaryChoice("att");
    S = new ActionCounter("sch",8,s);
    Actions(s);
    EndogenousStates(E);
    ⁞
    }

The user must tell niqlow to impose the condition that the agent can either work or study but not both. This creates an additional trimming of unreachable states: $M+S\le t.$ In this approach the agent has to impose this extra condition on reachable states. The user replaces two built-in virtual methods with their versions:

LSext::FeasibleActions() {
    return CV(m) .* CV(s) .== 0;
    }
LSext::Reachable(){
    return CV(M) + CV(S) <= I::t;
    }

Augmented State Variables

Here are 3 augmented state variables:

    x = new Freeze( new ActionCounter("sch",8,s), 15 );
    y = new Reset(b,a);
    z = new ForgetAtT( new ValueTriggered(d,tvar,1,5), 20 );

State variable $y$ augments a base state variable $b$ (not shown here) so that its value resets to 0 whenever the agent sets the action variable $a$ to 1 (also defined elsewhere). This is a special case of the general Triggered() augmentation. Several triggers besides the simple reset are already coded. Finally, z is a double augmentation of a state variable $d$ and another state variable $tvar.$ . First, when tvar=1 $z$ will reset to 5. From $t=20$ and onward the value of $z$ is not tracked because of the ForgetAtT augmentation. Forgetting a sate variable means its value is simply $CV(z)=0$ from then on, which avoids expanding the state space unnecessarily.

Complete Code and Output for the Labor Supply Example

#import "niqlow"
class LS : ExtremeValue {
    static decl m, M, e, beta, pi;
                Utility();
    static      Build(d=FALSE);
    static      Create();
    static      Earn();
    }
main() {
    LS::Create();
    VISolve();
    ComputePredictions();
    }
LS::Build(d) {
     SetClock(NormalAging,40);
     if (isint(d)) {
        e = new Nvariable ("e",15);
        m = new BinaryChoice("m");
        Actions(m);
        ExogenousStates(e);
        }
     else
        m = d;
     M = new ActionCounter("M",40,m);
     EndogenousStates(M);
     SetDelta(0.95);
     beta =<1.2 ; 0.09 ; -0.1 ; 0.2>;
     pi = 2;
    }
LS::Create() {
    Initialize(1.0,new LS());
    Build();
    CreateSpaces();
    }
LS::Earn() {
    return  exp( (1~CV(M)~sqr(CV(M))~AV(e)) * CV(beta) ) ;
    }
LS::Utility() {
    return CV(m)*(Earn()-pi) + pi;
    }

 t           m           M                 t           m           M
  0      0.3982      0.0000                20      0.1540      4.9016
  1      0.3761      0.3982                21      0.1508      5.0556
  2      0.3537      0.7743                22      0.1480      5.2064
  3      0.3319      1.1279                23      0.1456      5.3545
  4      0.3111      1.4598                24      0.1434      5.5000
  5      0.2918      1.7709                25      0.1414      5.6434
  6      0.2741      2.0628                26      0.1397      5.7848
  7      0.2580      2.3369                27      0.1382      5.9246
  8      0.2434      2.5949                28      0.1369      6.0628
  9      0.2303      2.8383                29      0.1358      6.1997
 10      0.2186      3.0686                30      0.1348      6.3355
 11      0.2082      3.2872                31      0.1339      6.4703
 12      0.1989      3.4954                32      0.1332      6.6042
 13      0.1907      3.6944                33      0.1325      6.7373
 14      0.1834      3.8850                34      0.1320      6.8699
 15      0.1769      4.0684                35      0.1316      7.0019
 16      0.1712      4.2453                36      0.1313      7.1335
 17      0.1661      4.4165                37      0.1310      7.2648
 18      0.1615      4.5825                38      0.1309      7.3958
 19      0.1575      4.7441                39      0.1308      7.5267
 

In this model we can push the smoothing parameter from 1.0 to near 0 and to very large:

    Initialize(1.0,new LS());         Baseline
    Initialize(0.001,new LS());       Near perfect smoothing
    Initialize(100.0,new LS());       Almost no smoothing

    SetDelta (0.95);                  Baseline
    SetDelta (0.0);                   Static decisions

    beta = <1.2 ; 0.09 ; -0.1 ; 0.2>;   Baseline
    beta = <1.2 ; 0.0; 0.0; 0.2>;       Static environment (no experience)

Unlike purpose-built code, in niqlow these kinds of tests use the same underlining code for all models. In addition, the behavior of a state variable class, such as ActionCounter can be confirmed in a small test program. Then it is very likely it will perform correctly in any other problem.

Certainly not all bugs have been discovered let alone fixed in the current niqlow code. And varying parameters of a single model will not reveal all errors. However, as changes are made that might break code that used to work, a suite of test programs are run to make sure that output is still correct. In that sense the changes to experiments on the labor supply model are important for a user to run in order to overcome healthy skepticism. However, the kinds of errors that might uncover have been squashed by checking test program output that include many other features than the labor supply model but in smaller spaces. There are also debugging features that can be turned on to trace output when a bug has been discovered.

Appendix C: Solution Algorithms

Appendix C: Reservation Values

Once $z^\star$ has been found, Bellman iteration requires calculation of the expected value of arriving at the (now explicit) state $\th:$ $$\eqalign{ EV(\th) = G(z^\star)&\l( E\l[U|a=0,z\lt z^\star\r] + \delta EV_0 \r)\cr + \l(1-G(z^\star)\r)&\l( E\l[U|a=1,z\ge z^\star\r] + \delta EV_1 \r).\cr\label{EVz} }$$ These conditions are illustrated in A2. Conditions \eqref{WS1}-\eqref{EVz} include three additional objects that the user's code must provide. The first condition for $z^\star$ needs $U\l(A(\theta),z\r),$ a vector of utilities for candidate values of $z$ while solving for $z^\star.$ The second includes $F(z^\star)$ and the vector of expected conditional utilities, $E\l[U\l(0\r)| z\lt z^\star\r]$ and $E\l[U\l(1\r)|z\ge z^\star\r].$

The model must satisfy several conditions described there to be eligible for this solution algorithms. These are enforced by allowing only models derived from the OneDimensionalChoice class to use the ReservationValues method.

1. $\al$ must be one-dimensional (only one variable added to the action vector);
2. No smoothing shock $\zeta$ is included (because $Z$ plays that role);
3. No state variables are placed in $\eps$ and $\eta$ vectors (because reservation values must be stored using $\th$);
4. Choice values $\vv$ must the single-crossing property.

OneDimentionsalChoice adds two virtual functions that other Bellman classes do not. The two functions are Uz() and EUtility(). At $z^\star$ the difference in current utility equals the discounted difference in future values. The user codes Uz(z) to return utility of all choices at $z$. The reservation value solution method uses that to solve for $z^\star.$ Backward induction requires the expected utility of each option conditional given the optimal value $z^\star$, hence the need to provide EUtility().

---

---

Code for the Reservation Value Labor Supply Model

Since $E$ exceeds the value of not working ($\pi$) for large enough values of $z$ there is a reservation value for the choice $m=1.$ Since LSz must be a one-dimensional choice model it cannot be derived from LS as was done earlier with LSext.

However, code from the basic model can be reused because most elements of LS are static. The only non-static member was Utility(). By using the static elements a change to LS will still be reflected automatically in LSz. In particular, LSz::Uz() can use LS::Earn() to utility at $z$. First set the value of LS::e to z then calls LS::Earn().

In the model expected utility of not working is simply $\pi$. For working, expected utility involves $E[e^{\beta_3z}|z>z^\star]$. The Mills-ratio formula for log-normality gives expected utility conditional of acceptance as $$E[U|m=1,z>z^\star] = \exp\l\{\beta_0+\beta_1M+\beta_2M^2+{\beta_3^2\over 2}\r\}{\Phi\l(z^\star/\beta_3 - \beta_3\r)\over \Phi\l(z^\star\r)}.\label{EU}$$

The constant factor includes $\beta_3^2/2$, because it is coefficient on the normal random variable $z$ (and hence the variance of the model shock). The ratio of $\Phi()$ values is new to the continuous specification and can't be borrowed from the discretized LS.

However, note that only the constant term includes the Mincer equation, and it is the same as the original earnings function if $e=\beta_3/2.$ So again the base LS::Earn() function can be used by setting the value of e first. Thus, even though LS was coded for a completely different approach its specification is still synchronized with the reservation wage version. Only elements specific to the new version need to be coded.

This code segment converts the labor supply model to a continuous choice reservation value problem.

#include "LS.ox"

class LSz : OneDimensionalChoice {
    static Run();
    EUtility();
    Uz(z);
    }
LSz::Run() {
    Initialize(new LSz());
    LS::Build(d);
    CreateSpaces();
    RVSolve();
    ComputePredictions();
    }
LSz::Uz(z) {
    LS::e = z;
    return LS::pi | LS::Earn();	
    }
LSz::EUtility()    {
    decl pstar = 1-probn(zstar),
         sig   = LS::beta[3];
    LS::e = sig/2;
	return {  ( LS::pi | LS::Earn()*probn((zstar/sig-sig)/pstar))
              , (1-pstar)~pstar
            };
	}	

Appendix D: Estimation on Simulated Labor Supply Data

    dta = new OutcomeDataSet("data",vi);
    dta -> Simulate(1000,40);
    dta -> ObservedWithLabel(m,M,obsearn);

Abbreviated output is below. The code detects that none of the observations are full CCP because $e$ is unobserved. Since all actions and endogenous states $\th$ are observed all paths are categorized as IID. This means niqlow will automatically sum over $\eps$ to compute the likelihood but it is unnecessary to use the backward Algorithm 6. The measurement error on noisy earnings is fixed at its true value. The result is convergence after 5 BHHH iterations with weak convergence. Because by default niqlow scales starting parameters the reported standard errors would need to be re-scaled, or recomputed with scaling and constraining turned off (details available in the documentation). Having started at the true parameter values, the sample likelihood is within .01% of the initial value with some differences in parameter estimates, particularly the coefficient on experience ($M$).

Masking data for observability.
Path like type counts
    CCP    IIDPartObs
      0   1000      0
Report of Gradient Starting on lnklk
   Obj=          -29378.9180146
Free Parameters
         index          free
B0         0     1.00000000000
B1         1     1.00000000000
B2         2     1.00000000000
B3         3     1.00000000000
Actual Parameters
                  Value
B0         1.20000000000
B1       0.0900000000000
B2       -0.100000000000
B3        0.200000000000

1. f=-29375.8 deltaX: 0.32545 deltaG: 50.0006
        Output for iterations 1-4 removed 
5. f=-29375.8 deltaX: 4.45343e-05 deltaG: 0.00906749

Finished: 3:WEAK
                          B0           B1           B2           B3
    Free Vector       1.0485      0.66901      0.97388      0.93059
    Gradient     -3.3945e-05  -4.1913e-05   0.00017144   1.0543e-05
   Std.Error        0.022743      0.14097     0.014054      0.18729

Report of  Iteration Done  on lnklk

   Obj=          -29375.7501799
Free Parameters
         index          free                  stderr
B0         0     1.04850905614   0.0227425353445
B1         1    0.669005636205    0.140973390182
B2         2    0.973884754314   0.0140539379652
B3         3    0.930592373935    0.187285663731
Actual Parameters
                  Value
B0         1.25821086737
B1       0.0602105072584
B2      -0.0973884754314
B3        0.186118474787

Declaration

Funding and/or Conflicts of interests

References