How can I estimate the parameters of an option value model of retirement?

I am modelling an option value model of retirement, see for instance Stock and Wise (1990). I am however not sure to which class of problems this model falls into and hence which optimization method I should consider to solve for the issue. To clarify what I mean, I will first describe the model:

The objective is to find probabilities of retiring in a specific year, and the parameters $\rho$, $\gamma$, $k$, $\beta$ in: \begin{eqnarray} Pr[\text{retire}_t] & = & Pr[ g_t( r^{*}_t ) / K_t( r^{*}_t ) < - v_t ], \\ v_{s} & = & \rho v_{s-1} + \epsilon_t \\ g_t(r_t) & = & \sum^{r-1}_{s=t} \beta^{s-t} \pi(s|t) E_t(Y_s^{\gamma}) + \sum^{S}_{s=r} \beta^{s-t} \pi(s|t) E_t( [k B_s(r)]^{\gamma} ) \\ & - & \sum^{S}_{s=r} \beta^{s-t} \pi(s|t) E_t( [ k B_s(t)]^{\gamma} ), \\ K_t( r_t ) & = & \sum^{r-1}_{s=t} \beta^{s-t} \pi(s|t) \rho^{s-t}, \\ v_t & = & (\omega_t - \xi_t ) \end{eqnarray} Here $Y_s$ are future wages and $B_s(t)$ are retirement incomes with $\pi(s|t)$ the probability that a person will live in year s given that she or he lives in year $t$. $r^{*}_t$ is the year in which the value of future stream of income is maximized. The value of future stream of income if retirement is at age $r$ is given by: \begin{eqnarray} V_t(r) & = & \sum^{r-t}_{s=t} \beta^{s-t} U_w(Y_s) + \sum^{S}_{s=r} \beta^{s-t} U_r[B_s(r)], \\ U_w(Y_s) & = & Y_s^{\gamma} + \omega_s, \\ U_r(B_s) & = & [ k B_s(r)]^{\gamma} + \xi_s, \\ \omega_s & = & \rho \omega_{s-1} + \epsilon_{\omega,s} \\ \xi_s & = & \rho \xi_{s-1} + \epsilon_{\xi, s} \end{eqnarray}

Usually one could solve such issues by considering it as a dynamic factor model or with maximum likelihood. However the summation $\sum$ makes it tricky as the 'optimal' date of retirment $r$ is a parameter itself.

So in short, how would you classify this problem and which optimization techniques do you recommend? Do you perhaps know packages or know where I can find a code to solve for this?

-

Have you tried a swaption having a marked to market of the notional N at t_1=retirement age?

From t_0=contract signing untill t_1, you receive the rate r_payments and pay the risk-free/investment rate. Then you pay out a predifined rate r_retirement and you pay in the risk-free/investment rate untill t_2=death time.

From your point of view, you are exposed to the default risk of the counterparty (death) given by the survival rate, with the recovery notional being N. If you do not invest at the risk-free rate, you are sensitive, in addition, to the credit risk of the counterparty you invest in. You are generating the credit risk of (defaulting) yourself.

All the entering parameters are practically known, especially if you do not make use of the risk-free interest rate, with credit risk modellable via the spreads.

The notional N could, in addition be a nicer function which, due to missing payments, is not necessarily proportional with (t_1-t_0) and which you might overpass also by an adjustment to r_payments.

-

In principle, this seems to be a job for policy iteration. In options pricing, we can use policy iteration for early exercise conditions (or other intermediate stopping times) as seen, for example, here and here.

In practice, since your problem of identifying $r^*$ is just one-dimensional and integer-valued, you may find it quicker and easier to simply iterate potential values of $r^*$ to identify the optimal one as an inner loop of your optimization process for $\rho, \gamma, k, \beta$.

This is to say, the code for the objective function you fit in determining the parameters looks like this:

myObj(rho, gamma, k, beta, observedProb) {
for r in (1:100) {
retireValue[r] = V_t(r)
}
rStar = argmax( retireValue )
dist = modelRetirePro(rStar, rho, gamma, k, beta) - observedProb
return dist*dist
}

-