Past exam of the mathematics course of the University of Cambridge / 2008 / iii / Paper 47 (Time Series and Monte Carlo Inference)

A weakly stationary process has finite second moments, a mean $m=E{X}_t$ independent of $t$, and a covariance depending only on the separation of the observations. Its autocovariance function and autocorrelation function arewhere the autocorrelation requires $\gamma (0)>0$.

Write the autoregressive polynomial as $A(z)=1-\alpha z+\alpha z^2$, so the autoregressive model is $A(B)X_t={ϵ}_t$, with $B$ the backshift operator. It is important to distinguish a weakly stationary process from a causal time series: weak stationarity alone does not require $X_t$ to depend only on present and past white noise.

For a two-sided weakly stationary process, the two-sided stationary inverse of an autoregressive polynomial exists and is unique if $A$ has no root on the unit circle. To find the excluded parameters, set $z=e^{i\lambda }$. The imaginary part of $A(z)=0$ givesThe case $\alpha =0$ is harmless. The possible frequencies are $0,\pi ,\pm \pi /3$. Since $A(1)=1$, $A(-1)=1+2\alpha$ and $A(e^{\pm i\pi /3})=1-\alpha$, the exceptional parameters are $-1/2$ and $1$. Thus, with nondegenerate white noise, the literal two-sided answer isFor completeness, away from those values $1/A(z)$ has an absolutely summable Laurent series on an annulus containing the unit circle. Its bilateral coefficients define an $L^2$-convergent linear filter of a stationary time series applied to ${ϵ}_t$, giving the required solution. Applying the same inverse to the equation gives uniqueness. At an exceptional value, the spectral measure of a stationary time series would have to satisfy $|A(e^{i\lambda }){|}^2{\mu }_X(d\lambda )={\sigma }^{2}d\lambda /(2\pi )$. A zero of $A$ on the unit circle makes the resulting density nonintegrable near that zero, contradicting finite variance.

If the usual additional causal time series convention is intended, the causality root criterion for an autoregressive model requires the roots of $A$ to lie outside the closed unit disk. Equivalently, the roots of $u^2-\alpha u+\alpha$ lie strictly inside it. The real quadratic stability conditions reduce to $|\alpha |<1$ and $1+2\alpha >0$, giving the narrower answerThe printed question does not explicitly impose this extra convention. Both answers have therefore been distinguished. If ${\sigma }^2=0$, the exceptional parameters permit nonunique stationary homogeneous solutions rather than the nonexistence conclusion above.

At $\alpha =-1/12$,Consequently the Wold representation isThe coefficients are absolutely summable, and ${\psi }_0=1$. The past of $X$ lies in the closed linear span of the past white noise; conversely ${ϵ}_t=A(B)X_t$ recovers that noise from present and past $X$. Thus ${ϵ}_t$ is orthogonal to the past of $X$ and is its one-step linear prediction error, as required for the Wold decomposition. There is no deterministic component.

Using orthogonality of the white noise, put $H=|h|$ and calculate the autocovariance by geometric sums:HenceIn particular, $\gamma (0)=66{\sigma }^2/65$ and $\gamma (1)=-6{\sigma }^2/65$; these also satisfy the Yule-Walker equations for this autoregressive model.

Use the two-sided spectral density of a stationary process conventionThus white noise of variance ${\sigma }^2$ has density ${\sigma }^2/(2\pi )$.

Absolute summability makes the linear filter of a stationary time series well defined in $L^2$: ${\sum }_r\Vert a_r{X}_{t-r}{\Vert }_2\le \Vert X_0{\Vert }_2{\sum }_r|a_r|<\infty$. Its mean is $m_Y=m_X{\sum }_r{a}_r$. For its autocovariance, continuity of the $L^2$ inner product givesThe sum is absolutely convergent because $|{\gamma }_X(j)|\le {\gamma }_X(0)$. It depends only on $k$, proving weak stationarity. Inserting the spectral inversion formula and exchanging the absolutely dominated sum with the integral yieldsHere the coefficients are real, so the two factors are complex conjugates. This is the spectral density transformation under a linear filter.

Applying a second linear filter of a stationary time series givesThe composition of absolutely summable time-series filters can also be calculated directly:Indeed ${\sum }_{\ell }|c_{\ell }|\le \Vert a{\Vert }_1\Vert b{\Vert }_1$; this bound also justifies the interchange of the two $L^2$ sums. The filter gain multiplies under this composition. The sketches for the two individual filters are:

For differencing, the only nonzero filter coefficients are $a_0=1$ and $a_1=-1$. Therefore $a(\lambda )=1-e^{i\lambda }$ andThe left sketch shows the filter gain increasing from zero at frequency zero to two at $\pi$. This differencing filter removes a constant component and attenuates slowly varying components while relatively emphasizing higher frequencies. It is a high-pass filter, rather than an averaging or smoothing filter.

The seasonal difference operator has $b_0=1$, $b_{12}=-1$ and all other coefficients zero. Its frequency response is $1-e^{12i\lambda }$, soThe right sketch shows zeros at $\lambda =j\pi /6$ for $j=0,\dots ,6$, and maxima of height two at $(2j+1)\pi /12$ for $j=0,\dots ,5$. The seasonal difference operator removes every exactly twelve-periodic component, including a constant. Its filter gain has a comb pattern: it selectively suppresses seasonal frequencies, rather than increasing monotonically with frequency.

The composition of absolutely summable time-series filters is $(1-B^{12})(1-B)$, with $B$ the backshift operator. Thusand its only nonzero coefficients are $c_0=1$, $c_1=-1$, $c_{12}=-1$, $c_{13}=1$. Multiplication of the two filter gains givesEquivalently, its spectral density transformation under a linear filter is $f_Z=16{\sin }^2(\lambda /2){\sin }^2(6\lambda )f_X$.

A good pseudorandom number generator for Monte Carlo integration should have approximately uniform outputs, a period much longer than the intended run, and negligible detectable dependence both serially and in higher-dimensional tuples. It should be fast, reproducible from a seed, portable enough to reproduce a calculation, and numerically precise enough for the required tail probabilities. Values strictly inside $(0,1)$ avoid singularities in logarithmic inverse transform sampling. Separate simulations also need well-separated streams or seeds, rather than accidentally identical streams. A pseudorandom number generator is deterministic; these properties mean that it imitates independent uniform distributions sufficiently well for the application, not that its outputs are literally independent.

For rate $\lambda >0$, the exponential distribution has cumulative distribution function $F(x)=1-e^{-\lambda x}$ on $x\ge 0$. Inverse transform sampling gives $-\log (1-U)/\lambda$; since $1-U$ is also uniform, an equivalent answer isFor $x\ge 0$, $P(E>x)=P(U<e^{-\lambda x})=e^{-\lambda x}$, verifying the exponential distribution directly.

The centered Laplace distribution with density $(\lambda /2)e^{-\lambda |x|}$ can be obtained by multiplying $E$ by an independent sign taking $1$ and $-1$ with equal probabilities. Alternatively, its inverse transform sampling formula uses just one uniform variate:The first branch inverts the negative half of the cumulative distribution function and the second branch inverts the positive half. Add a location parameter if a noncentered Laplace distribution is wanted.

Correct uniformity is essential for each marginal distribution, and independence is essential for independent samples and for the sign construction. For example, setting every $U_i$ equal to a single uniform $U$ gives uniform marginals but identical transformed draws. Likewise, choosing a sign from the same uniform used for the magnitude can change the Laplace distribution unless the one-uniform inverse above is used.

The first algorithm uses independent normal distributions: generate $Z_1,\dots ,Z_{\nu }$ with law $N(0,1)$ and returnFor example, the Box-Muller transform generates two such normal variates from independent uniforms:To verify this, $R^2=-2\log U$ has density $e^{-s/2}/2$ and $2\pi V$ is an independent uniform angle. Converting the resulting polar density to Cartesian coordinates gives $(2\pi {)}^{-1}{e}^{-(z_1^2+z_2^2)/2}$, which factors into two independent normal distributions. Independent pairs therefore give the required sum, discarding a spare coordinate when $\nu$ is odd. Its moment-generating function is $(1-2t{)}^{-\nu /2}$, that of the chi-squared distribution.

A genuinely different algorithm uses exponential-envelope rejection sampling for a chi-squared variable. Put $a=\nu /2\ge 1$. The target and proposal probability density functions areFor $a>1$, differentiating $\log (f/g)$ gives $(a-1)(1/y-1/\nu )$, so the maximum is at $y=\nu$. The envelope constant and normalized density ratio areThey remain valid at $a=1$, where $M=1$ and $f=g$. Repeatedly draw fresh independent $U,V$, put $Y=-\nu \log U$, and accept whenThe right side is nonpositive by $\log u\le u-1$, so this is a legitimate acceptance probability. The joint probability of proposing and accepting in $dy$ is $g(y)f(y)/(Mg(y))dy=f(y)dy/M$; conditional on acceptance the density is exactly $f$. This proves the rejection sampling algorithm, rather than only identifying the target law.

For $\nu =2$ I prefer the second method: it reduces to a single exponential draw with no rejection. The first method is simple and attractive when normal draws are already available or $\nu$ is small. For larger $\nu$, the second method avoids generating $\nu$ normals: its acceptance probability is $1/M$, asymptotic to $\sqrt{2\pi }/(e\sqrt{a})$ by Stirling's formula. Its expected number of proposals grows like $\sqrt{\nu }$, rather than the linear number of normal draws in the first method. Both are exact ideal-uniform algorithms; actual runtime also depends on the implementation.

For Beta sampling by gamma ratios, use either preceding chi-squared distribution algorithm to generate independent $S\sim {\chi }_{2\alpha }^2$ and $T\sim {\chi }_{2\beta }^2$, and returnWrite $W=S+T$, so $S=WR$ and $T=W(1-R)$. The Jacobian determinant is $W$. The joint probability density function after this change of variables factors asIntegrating out $w$ proves the Beta distribution for $R$.

A new method is Beta sampling by uniform order statistics. Generate $m=\alpha +\beta -1$ independent uniforms, sort them, and return the $\alpha$th order statistic. There are $\alpha -1$ uniforms below the selected value and $\beta -1$ above it. Selecting their labels and placing the selected observation in $dr$ gives densitywhich is the Beta distribution density. The integrality of the shape parameters is what makes this order-statistic construction possible.

For jackknife resampling, calculate the leave-one-out estimates ${\hat{\theta }}_{(-i)}$ by deleting observation $i$, and their average ${\overline{\theta }}_J=n^{-1}{\sum }_i{\hat{\theta }}_{(-i)}$. The jackknife variance estimator isThe factor $n-1$ compensates for the small differences between the overlapping samples. For sufficiently smooth statistics it estimates their sampling variance. The related estimated bias is $(n-1)({\overline{\theta }}_J-\hat{\theta })$, leading to the jackknife bias correction $n\hat{\theta }-(n-1){\overline{\theta }}_J$.

For the nonparametric bootstrap, replace the unknown distribution by its empirical distribution functionIndependently sample $n$ observations with replacement from the original observations, compute the estimator on that bootstrap sample, and repeat this for $B$ independently generated bootstrap samples. If the resulting values are ${\hat{\theta }}_1^{\ast },\dots ,{\hat{\theta }}_B^{\ast }$, estimate the conditional bootstrap variance byIts interpretation as an approximation to the original sampling variance relies on the bootstrap approximating the sampling law of the statistic; increasing $B$ only reduces simulation error and cannot repair a failure of that approximation.

Since the original observations are distinct, a bootstrap sample has no repeated value precisely when it contains all $n$ original values once. Of the $n^n$ equally likely ordered resamples, $n!$ are permutations without repetition. ThereforeThe probability is zero for $n=1$ and rapidly approaches one as $n$ grows.

The first matrix contains the leave-one-out samples, and each element of vect is their sample mean. For $n\ge 2$,Thus R4a computes the jackknife variance estimator of the full sample mean:where $s_x^2$ is the unbiased sample variance. R5a computes the estimated bias used in jackknife bias correction. The identity for the average leave-one-out sample mean gives R5a = 0, apart from floating-point rounding.

The second printed loop contains an indexing error: its row index is i, although its loop variable is b. Literally, after the first block, i equals $n$. When $n\le 199$, every iteration overwrites row $n$ of the new matrix; the other rows remain missing. Hence vect contains one finite mean and 198 missing entries. R's default sorting removes missing entries, leaving a vector of length one, so selecting its fifth and 195th entries gives R7b = c(NA, NA). If $n>199$, the matrix assignment is out of bounds; if the second block runs independently with i undefined, it fails even earlier. This is a source error, not a property of the bootstrap.

Replacing just the row index by b gives the intended algorithm: draw 199 independent bootstrap samples and sort their sample means as $s_1\le \cdots \le s_{199}$. The two indices are $200(0.05/2)=5$ and $200(1-0.05/2)=195$, soThese are the endpoints of an approximate 95% percentile bootstrap confidence interval for the population mean. Their actual numerical values depend on the supplied data and random draws, which are not specified. They are not the reflected endpoints of a basic bootstrap confidence interval, and the choice of ranks does not guarantee exact finite-sample coverage.

Negative covariance can be exploited through antithetic variates, but it must belong to an actual joint sampling scheme. Two separately independent samples do not inherit negative covariance merely because one can imagine a coupling of their marginal distributions. Interpret the data as $h=\min (n,m)$ independent matched pairs $(X_i,Y_i)$ with covariance ${\sigma }_{XY}<0$, together with independent unmatched observations if the sample sizes differ. Other known pairing schemes can be handled by substituting their actual covariance of the two means.

Let $A=\overline{X}$, $B=\overline{Y}$, and putFor any fixed $w$, the combination $wA+(1-w)B$ is an unbiased estimator of the common mean. Its variance is $w^2{V}_X+(1-w{)}^2{V}_Y+2w(1-w)C$. Differentiation, or completing the square, gives the minimum-variance combination of unbiased sample means:With $C<0$, both weights are positive. This is efficiency among fixed linear unbiased combinations, not a universal claim about every possible estimator for arbitrary distributions $F,G$.

For an immediately implementable exactly unbiased estimator without known moments, use fixed $w=n/(n+m)$, the pooled sample mean. Its variance iswhich shows directly the gain from negative covariance. For any fixed weight, estimate $V_X,V_Y$ using $s_X^2/n,s_Y^2/m$, and $C$ using $h{s}_{XY}/(nm)$, with $s_{XY}$ the usual unbiased sample covariance of the matched pairs. Substitute them into the variance quadratic. At least two matched pairs are needed to estimate their covariance nonparametrically.

If the unknown optimal weights are desired while retaining exact unbiasedness, a cross-fitted combination of unbiased sample means gives an algorithm using the supplied data. Split the independent pairs and unmatched observations into two folds, keeping each pair together and retaining observations of both types in each fold. Estimate the moments from fold 2 to choose the weight for the fold-1 sample means, adjusting the variance terms for fold 1's sample sizes; do the converse for fold 2. Clip estimated weights to $(0\le w\le 1)$, to keep them bounded when training moment estimates are unstable. Average the two held-out combinations with deterministic weights, for example one half each for balanced folds.

Conditional on its training fold, each held-out pair of sample means has expectation $\theta$, and its weights sum to one. Each held-out combination, and hence their average, is therefore exactly an unbiased estimator. Under consistent moment estimation with increasing balanced folds, the weights approach $w_{\ast }$ and recover its first-order variance advantage. In contrast, estimating a weight from the same data to which it is applied need not yield an unbiased estimator.

Estimate the variance of this fitted-weight procedure by a bootstrap that resamples matched pairs as blocks and unmatched observations separately, repeats the entire split/weight/combination procedure, and takes the sample variance of the resulting estimates. This preserves the negative covariance and includes variability of the fitted weights. Resampling the $x$ and $y$ values independently would erase the feature responsible for the improvement. The bootstrap estimate has its usual large-sample qualification; the fixed-weight plug-in variance formula above is the simpler alternative when weights are predetermined.

A systematic Gibbs sampler updates the coordinates of $\theta$ successively from their full conditional distributions. At iteration $t$, drawFor a single coordinate, if the current state has distribution $\pi$, its untouched coordinates have the correct marginal law and the resampled coordinate has precisely the conditional law under $\pi$. Integrating marginal times conditional therefore recovers $\pi$. Each coordinate kernel preserves this invariant distribution, so their composition also preserves it. This proves the target of the Gibbs sampler; the whole systematic sweep need not be a reversible kernel.

Starting from a suitable state and assuming the required irreducibility and convergence conditions, discard an initial transient and retain subsequent states as a dependent sample from, or approximately from, $\pi$. In finite-state problems irreducibility and aperiodicity suffice; in continuous-state problems the corresponding positive recurrence and ergodicity conditions must hold. Successive draws are generally dependent, so Monte Carlo uncertainty must account for that dependence.

Use the shape-scale convention for the inverse-gamma distribution: $\mathrm{IG}(a,b)$ has density proportional to $s^{-a-1}{e}^{-b/s}$ for $s>0$. Write $X=X_k$, $\beta ={\beta }_k$, $\mu ={\mu }_k$, $\Sigma ={\Sigma }_k$ and $s={\sigma }^2$. The independent normal and inverse-gamma regression priors make the conditional log-likelihood plus log prior in $\beta$ equal, up to a constant, toDefine $V=(X^{T}X/s+{\Sigma }^{-1}{)}^{-1}$ and $m=V(X^{T}y/s+{\Sigma }^{-1}\mu )$. Completing the square turns this expression into $-(\beta -m{)}^T{V}^{-1}(\beta -m)/2$ plus a constant. Thus the conditional posterior distribution is the multivariate normal distributionPositive definiteness of the prior covariance matrix makes this well defined even if the design matrix is rank deficient.

Holding $\beta$ fixed, the factors depending on $s$ areConsequentlyThere is no additional $(k+1)/2$ in the shape: the normal distribution prior for ${\beta }_k$ is independent of ${\sigma }^2$. A prior scaled by ${\sigma }^2$ would be a different model and would change that calculation.

Choose an initial $s^{(0)}>0$. At each iteration, evaluate $V$ and $m$ using $s^{(t)}$ and draw ${\beta }^{(t+1)}$ from the conditional multivariate normal distribution just derived. Next compute $b_t=b+\Vert y-X_k{\beta }^{(t+1)}{\Vert }^2/2$ and drawFor example, generate $G\sim \Gamma (a+n/2,\text{rate }{b}_t)$ and set $s^{(t+1)}=1/G$. A draw of the conditional multivariate normal distribution can be formed as $m+LZ$ with $L{L}^T=V$ a Cholesky decomposition and $Z$ a vector of independent standard normal distributions.

This two-block Gibbs sampler alternates the full conditional distributions of the independent normal and inverse-gamma regression priors posterior. The conditional-integration argument proves that its invariant distribution is the joint posterior distribution. Under the usual convergence conditions, retain the pairs $({\beta }^{(t)},s^{(t)})$ after an initial transient as the requested dependent posterior sample. The new variance update uses the newly drawn coefficients, rather than the previous coefficients.

A reversible-jump Markov chain Monte Carlo calculation compares normalized model-specific densities. The fixed-order posterior written only up to proportionality in part (a) omits normalizing factors that depend on $k$ and therefore cannot simply be canceled between dimensions. Moreover, the printed question does not specify the model-order prior distribution or the probabilities of selecting birth and death moves. Denote those by $p_k$, $b_k$ and $d_k$, respectively, retaining them in the answer.

Put $s={\sigma }^2$, $c=n^{-1}{\sum }_i{x}_i^{k+1}$, $v_i=x_i^{k+1}-c$, and $r=y-X_k\beta$. The centered polynomial birth move in reversible-jump sampling changes the fitted vector by $zv$, so its new residual is $r^{\prime }=r-zv$. ThereforeThe map from $({\beta }_0,\dots ,{\beta }_k,z)$ to $({\beta }_0^{\prime },\dots ,{\beta }_{k+1}^{\prime })$ has ones down the diagonal and a single additional entry $-c$ in the intercept row. Its Jacobian determinant is one.

Write $Q_k=(\beta -{\mu }_k{)}^T{\Sigma }_k^{-1}(\beta -{\mu }_k)$ and $Q_{k+1}=({\beta }^{\prime }-{\mu }_{k+1}{)}^T{\Sigma }_{k+1}^{-1}({\beta }^{\prime }-{\mu }_{k+1})$. The normalized multivariate normal distribution prior ratio isDivide also by the forward proposal density $q(z)=(2\pi {\sigma }_{\beta }^2{)}^{-1/2}{e}^{-z^2/(2{\sigma }_{\beta }^2)}$. The common inverse-gamma distribution prior for $s$ and the common normal-likelihood normalization cancel, since $s$ stays fixed. The forward acceptance ratio is consequentlyHere ${\sigma }_{\beta }>0$ is the proposal standard deviation. If equal model priors and symmetric birth/death selection are explicitly adopted, their ratios become one; they are not determined by the PDF. The expression assumes the same hyperparameters for the variance prior in both models, as in part (a). A dimension-dependent variance prior would contribute its own ratio.

For the reverse centered polynomial birth move in reversible-jump sampling, recover $z={\beta }_{k+1}^{\prime }$ from the current coefficient vector. Remove that highest-order coefficient and undo the intercept shift:No auxiliary random number is drawn for this death move: the removed coefficient supplies the auxiliary value needed to reconstruct the birth move. The reverse Jacobian determinant is also one. The inverse prediction change is $-zv$, hence the residual in the lower-order model is $r=r^{\prime }+zv$.

Reversible-jump Markov chain Monte Carlo requires the reciprocal density/proposal ratio for the matched reverse transformation, so the death acceptance probability isExpressed directly in the current residual $r^{\prime }=y-X_{k+1}{\beta }^{\prime }$ and the recovered lower-order parameters, this isIndeed $r=r^{\prime }+zv$ makes this exponent the negative of the birth exponent, and all prefactors invert. Thus $R_{-}{R}_{+}=1$ for matched states. This reciprocity, including the normalized prior distributions and auxiliary proposal density, gives the required detailed balance between models.

The expectation-maximization algorithm alternates conditional averaging over missing data with maximization over the parameter. Start at a parameter value ${\theta }^{(0)}$ for which the observed likelihood function is positive. At iteration $t$, its E-step formsusing the current conditional law of $Z$ but treating $\theta$ as the candidate parameter. Its M-step takesRepeat until the parameter or observed log-likelihood stabilizes. An increase of $Q$ implies an increase of the observed log-likelihood: if $q_t(z)=f(z\mid x;{\theta }^{(t)})$, thenwhere Kullback-Leibler divergence is nonnegative. In general this monotonicity does not guarantee a global maximum; distinct initializations may approach different stationary points. The particular observed log-likelihood below is concave, which makes a stronger conclusion possible there.

The parameter-dependent part of the complete multinomial distribution log-likelihood isThe remaining probability factors and multinomial coefficient are constant in $\theta$. If $y_1+y_2+y_3>0$, differentiating gives $-y_1/(1-\theta )+(y_2+y_3)/\theta$, whose zero isThe log-likelihood is concave. If $y_2+y_3=0<y_1$, its maximum is at zero; if $y_1=0<y_2+y_3$, its maximum is at one. These boundary cases are included by the same formula. If $y_1=y_2=y_3=0$, the data lie entirely in the parameter-independent fourth cell and every $\theta \in [0,1]$ is a maximum-likelihood estimator.

Within the merged third observed cell, the probability of belonging to the original fourth cell isThe allocation therefore has the conditional binomial distributionIts conditional expectation and the filled-in missing counts arewith $y_1^{(t)}=x_1$ and $y_2^{(t)}=x_2$.

Taking the conditional expectation of the complete log-likelihood gives the EM for merged multinomial cells E-step:This is the parameter-dependent part of $\log L(y^{(t)}\mid \theta )$. More precisely, $Q(\theta \mid {\theta }^{(t)})=\log L(y^{(t)}\mid \theta )+C_t^{\prime }$, where the constant is independent of the candidate $\theta$, interpreting factorials of noninteger filled counts through the gamma function. The expected log-factorials are not generally the log-factorials of expected counts, so literal equality including the normalization would be false. Equality up to this constant is exactly what the M-step needs. Filling in conditional means works here because the candidate-parameter terms are linear in the missing counts; it is not a universal rule for the expectation-maximization algorithm.

Initialize $0<{\theta }^{(0)}<1$. At each iteration first compute $y_3^{(t)}=x_3{\theta }^{(t)}/(2+{\theta }^{(t)})$ by the conditional expectation above. The M-step applies the complete-data maximum-likelihood estimator to the filled counts, givingRepeat until convergence. If all counts are zero, the likelihood function is constant and any parameter can be retained. For a nonempty sample an interior start makes the displayed denominator positive.

Here the three observed cell probabilities are $(1-\theta )/2$, $\theta /4$ and $(2+\theta )/4$, so the observed log-likelihood isFor a nonempty sample its second derivative on $(0,1)$ isThus the maximum-likelihood estimator is unique, allowing a boundary maximum. To verify that the iteration actually finds it, let $u(\theta )=x_3\theta /(2+\theta )$ and $F(\theta )=(x_2+u(\theta ))/(x_1+x_2+u(\theta ))$. Then $F^{\prime }(\theta )=x_1{u}^{\prime }(\theta )/(x_1+x_2+u(\theta ){)}^2\ge 0$, andIf the maximum is interior, this identity moves each iterate towards it, and monotonicity of $F$ prevents crossing it. The bounded monotone iterates converge; continuity makes their limit a fixed point, and strict concavity makes that fixed point the maximum. If the maximum is at zero or one, the sign of ${\ell }^{\prime }$ gives the same monotone convergence towards the appropriate endpoint. In the special case $x_1=0$, $F$ is identically one whenever its denominator is positive, so the maximum is reached in one step. An interior initialization avoids the possible spurious zero fixed point when $x_2=0$; starting at zero can otherwise trap the iteration even when the true maximum is positive.