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

Interpret the forecasts as best linear prediction from a finite past. For a causal autoregressive process, the innovations are orthogonal to previous observations. If they also have zero conditional mean given the past, the same forecasts are conditional means; white noise by itself guarantees only the linear-prediction interpretation.

Let $m$ be the process mean and write its centered autoregression as $X_t-m={\sum }_{j=1}^p{\varphi }_j(X_{t-j}-m)+{ϵ}_t$. Set ${\hat{X}}_{T,h}=X_{T+h}$ for $h\le 0$, and recursively putThe forecast is a linear combination of the observed last $p$ values. Subtracting the recursion from the future autoregression shows that the prediction error is a linear combination of future innovations, hence orthogonal to every observed value. This proves the optimal projection property, rather than merely suggesting that future noise should be replaced by zero. In particular,The empty sum for $p=1$ is zero. Taking $m=0$ gives the centered convention. This is the recursive forecasts of an autoregressive process construction.

For the stable order-one equation, use its stationary causal solutionThe series converges in mean square because $\sum {\varphi }^{2j}<\infty$ and the white noise innovations have a common finite variance ${\sigma }^2$. Orthogonality of their distinct time indices givesAdding the deterministic mean therefore yields weak stationarity for $Z$, withThe autocorrelation function assumes ${\sigma }^2>0$; a zero-variance process has no normalized correlation. A recursion from an arbitrary nonstationary initial value would have transient moments, so the stationary-solution interpretation is important here.

Iteration gives $X_{T+k}={\varphi }^k{X}_T+{\sum }_{j=1}^k{\varphi }^{k-j}{ϵ}_{T+j}$. The innovation sum is orthogonal to the observed past, givingFor the integrated process, the backshift operator subtracts consecutive levels, so $(I-B)W_t=\mu +(I-B)Y_t=\mu +X_t=Z_t$. Consequently $W_{T+k}=W_T+{\sum }_{j=1}^k{Z}_{T+j}$. If the last increment is observed, substitute its forecast and sum the geometric progression:This is the forecasts of an integrated AR(1) process formula. Its final term is bounded as $k\to \infty$, even for negative $\varphi$, hence the forecast's long-run slope is the drift:With an initial level orthogonal to future innovations, the level formula is also the best linear forecast from the observed levels. Without that initial-level condition, it is the increment-based recursive forecast, and the levels could contain additional predictive information. The level formula uses $T\ge 2$, or a known $W_0$ when $T=1$. With just $W_1$ and no model for the initial level, the last increment cannot be recovered from the supplied observations; a one-observation level forecast would need that additional information.

For a weakly stationary process with mean $m$, the autocovariance function is ${\gamma }_k=\mathbb{E}[(X_{t+k}-m)(X_t-m)]$, independent of $t$. The positive definiteness of these covariances gives a unique finite positive spectral measure of a stationary time series on the frequency circle. Its cumulative function is the time-series spectral distribution $F$, with total mass ${\gamma }_0$, andUse one representative of the identified endpoints $-\pi$ and $\pi$ when specifying atoms. For a real-valued process the measure is symmetric, so the imaginary part integrates to zero. If it has a time-series spectral density, thenThe second identity is pointwise for an absolutely summable covariance sequence, as will apply to part (a). Existence of a density alone does not guarantee ordinary pointwise convergence of its Fourier series; in general the reconstruction can be interpreted by its Cesàro means in $L^1$. This fixes the frequency normalization for all three calculations below.

Only overlapping innovations contribute to the autocovariance function of the moving-average process. Set $c_0=1,c_1={\theta }_1,c_2={\theta }_2$. The overlap formula is ${\gamma }_k={\sigma }^2{\sum }_j{c}_j{c}_{j+|k|}$, with coefficients outside $0,1,2$ set to zero. ThusThe time-series spectral density is the Fourier sum of this finite covariance sequence:It also equals ${\sigma }^2|1+{\theta }_1{e}^{-i\lambda }+{\theta }_2{e}^{-2i\lambda }{|}^2/(2\pi )$, which directly verifies nonnegativity. Integrating from $-\pi$ gives the requested time-series spectral distribution:Extend it by zero below $-\pi$ and by ${\gamma }_0$ above $\pi$. In particular its total mass is the process variance, not necessarily one.

The even time-series spectral density makes the sine part vanish. At lag zero its integral is the area of the triangle, so ${\gamma }_0={\pi }^2$. For any nonzero integer $k$, integration by parts givesTherefore the triangular spectral density of a stationary time series has autocovariance functionNo additional factor $1/(2\pi )$ belongs in this inverse integral under the convention stated above.

The zero means of both amplitudes give $\mathbb{E}{X}_t=0$. Their unit variances and zero covariance give, for arbitrary integer $s,t$,Thus the variance is one and the covariance depends only on lag, proving weak stationarity. Independence or normality of the amplitudes is unnecessary. Sincethe spectral measure of a random harmonic oscillation assigns mass $1/2$ to each of $-{\omega }_0,{\omega }_0$. The corresponding right-continuous time-series spectral distribution isThese spectral atoms represent a persistent oscillation. There is no ordinary absolutely continuous spectral density for this process.

For rejection sampling, independently generate a proposal $Y$ with density $g$ and a uniform $U$ on $(0,1)$. Accept $Y$ if $U\le f(Y)/(Mg(Y))$; otherwise generate a fresh independent pair and repeat. Set the ratio to zero on $g=0$ where necessarily $f=0$. The envelope condition makes the acceptance threshold at most one.

For any measurable set $A$, a single attempt satisfiesTaking the whole space shows that acceptance has probability $1/M$. Conditioning on acceptance therefore gives density $f$. Repeating rejected independent attempts preserves that conditional law and terminates almost surely; the number of attempts is geometric with mean $M$.

For fixed proposal density, admissibility requires $M\ge {\sup }_{g(x)>0}f(x)/g(x)$ under the pointwise convention in the question. Since $1/M$ decreases with $M$, the smallest valid envelope gives the greatest acceptance probability:For densities specified only almost everywhere, replace the supremum by the essential supremum. Integration of the envelope also gives $M\ge 1$.

Use the uniform proposal on $[1,7]$, with density $g=1/6$. The triangular distribution peaks at $x=3$ with height $1/3$, so the optimal constant for this proposal is $M=2$. At each attempt take two fresh pseudo-random uniforms $U,V$, and propose $Y=1+6U$. The rejection method accepts according toIndeed these bounds are $f(Y)/(2g(Y))=3f(Y)$. Part (a) proves that the accepted value has the required density. Acceptance probability is $1/2$, so the method needs two proposal attempts on average, and hence four uniforms on average if each attempt uses two. The validity of the ideal algorithm presumes independent uniform draws; merely belonging to $[0,1]$ does not by itself make a pseudo-random sequence uniform or independent.

Integrating the two linear pieces gives the cumulative distribution functionThe split probability is $F(3)=1/3$. Invert each branch to obtain the method of inversion:Both expressions equal three at the split. For an ideal uniform draw, monotonicity gives $\mathbb{P}(F^{-1}(U)\le x)=\mathbb{P}(U\le F(x))=F(x)$, proving the target law. Inversion is preferable here: its quantile is explicit, it uses one uniform and one square root per output, and has no rejected draws. Rejection remains useful for densities with no convenient inverse CDF.

The nonparametric bootstrap replaces the unknown sampling law by the empirical distributionGenerate $B$ independent bootstrap samples, each containing $n$ independent draws with replacement from that law. Recompute the original estimator in each sample, obtaining ${\hat{\theta }}_1^{\ast },\dots ,{\hat{\theta }}_B^{\ast }$. With ${\overline{\theta }}^{\ast }=B^{-1}{\sum }_b{\hat{\theta }}_b^{\ast }$, the bootstrap standard error estimate isThis estimates the conditional bootstrap standard deviation, which approximates the estimator's sampling standard deviation when the bootstrap is consistent. Dividing this expression by $\sqrt{B}$ would instead estimate the Monte Carlo error in the bootstrap mean, a different quantity.

For distinct original observations, each resample is an ordered index sequence with one of $n^n$ equally likely values. It contains no repeats precisely when it is a permutation of all $n$ original indices, giving $n!$ possibilities. HenceThe distinctness assumption is needed to equate repeated values with repeated indices.

Calculate the observed sample correlation coefficient $\hat{r}$. In each paired bootstrap, draw 100 indices independently with replacement from $1,\dots ,100$ and select the corresponding whole pairs. Compute the bootstrap correlation $r_b^{\ast }$ from each resample; do not resample the two margins independently, since that would destroy their empirical dependence.

Let $q_p$ be the empirical $p$-quantile of those bootstrap correlations. The percentile bootstrap confidence interval isThis is the requested algorithmic interval; the actual endpoints require the observed pairs. Its coverage is approximate, based on the sampling law of the correlation being well approximated by the bootstrap law. Both sample variances must be positive. Resamples having a constant margin have undefined correlation and cannot be repaired by arbitrarily assigning correlation zero. One may regenerate these exceptional resamples, explicitly conditioning on well-defined correlations; under the usual nondegenerate large-sample conditions their probability is negligible. Frequent degeneracy instead indicates that this ordinary correlation-bootstrap construction is unsuitable.

For a bootstrap-t confidence interval, estimate the original standard error $\hat{s}$ of $\hat{r}$, for example by the standard deviation of the outer paired bootstrap correlations. For each outer resample $b$, estimate its own standard error $s_b^{\ast }$ by drawing many inner paired resamples from that outer sample and taking the standard deviation of their correlations. The resulting studentized statistics areLet $t_p^{\ast }$ be their empirical quantiles. Approximating the distribution of $(\hat{r}-r)/\hat{s}$ by that of $T^{\ast }$ and solving the two inequalities for $r$ givesThe reversal of quantiles is essential. The pivot requires positive standard errors and nondegenerate correlations; its claimed 95% coverage is an approximation under bootstrap regularity, not an exact finite-sample guarantee. One may intersect the final interval with the parameter space $[-1,1]$.

The same percentile and studentized procedures can be applied to the sample mean of $X$, and then only the $X$ observations need be resampled. They are useful when a normal approximation is doubtful. However the mean already has the simple standard-error estimate $S_X/10$ here. For normal observations, the classical exact interval isand the same form is an asymptotic approximation for independent finite-variance observations under appropriate large-sample conditions. Thus bootstrap intervals are valid alternatives, but are often unnecessary for a mean with a reliable classical pivot; the correlation has no corresponding distribution-free elementary pivot. An infinite-variance distribution would require other theory rather than this ordinary standard-error argument.

A control variate is a random variable $C$ correlated with the estimator and with known mean $c_0$. For a fixed real coefficient $a$, define $Y_a=Y-a(C-c_0)$. Then $\mathbb{E}{Y}_a=\theta$, whileAssuming finite second moments and positive control variance, complete the square to findWhen $\mathrm{Var}(Y)>0$, the minimum is $\mathrm{Var}(Y)(1-{\rho }_{Y,C}^2)$. A constant control cannot reduce variance. A strongly correlated control, including a negatively correlated one with a negative coefficient, can substantially reduce it.

For multivariate control variates, write $Y_{\beta }=Y-{\beta }^T(C-\mathbb{E}C)$, $\Sigma =\mathrm{Cov}(C)$ and $c=\mathrm{Cov}(C,Y)$. Completing the multivariate square givesIf $\Sigma$ is singular, use its Moore-Penrose inverse. Indeed any null direction has a constant centered control combination and therefore zero covariance with $Y$, so $c$ lies in the range of $\Sigma$. The elementary unbiasedness statement assumes a deterministic coefficient; estimating it from the same simulation may introduce finite-sample bias, whereas an independently trained coefficient keeps the conditional unbiasedness argument valid.

For the Metropolis–Hastings algorithm, start at a state $\theta$ with positive target density, propose ${\theta }^{\prime }$ from $q({\theta }^{\prime }\mid \theta )$, and generate an independent uniform draw. Accept the proposal with probabilityOn rejection retain $\theta$, including repeated states in the chain; removing repeats would generally change its stationary law. The accepted-flow density satisfieswhich is symmetric in its two states. This proves detailed balance and hence invariance of the target distribution, including the holding probabilities on the diagonal. Unknown common normalization constants cancel in the ratio. Invariance alone does not establish convergence from arbitrary initial states. Choosing proposals that yield a positive Harris recurrent Markov chain with aperiodicity supplies the usual convergence guarantee, and the resulting sample is dependent. Use burn-in and account for autocorrelation when assessing Monte Carlo precision.

A random-walk update proposes ${\theta }^{\prime }=\theta +\eta$, with an increment distribution independent of the current state. If its density is symmetric, $q({\theta }^{\prime }\mid \theta )=q(\theta \mid {\theta }^{\prime })$, reducing the Metropolis–Hastings acceptance probability to $\min \{1,\pi ({\theta }^{\prime })/\pi (\theta )\}$. For nonsymmetric increments the proposal-density ratio must remain.

Such updates need little global knowledge of the target and adapt naturally to its local scale. Their limitation is a tuning tradeoff: very small increments accept often but move slowly, while very large increments are often rejected. They can also have difficulty moving between distant modes. Covariance-scaled increments help with anisotropic targets, but the chain still explores through local steps.

An independence sampler proposes from a fixed density $g({\theta }^{\prime })$, independent of the current state. Its acceptance probability isIf $g$ is close to the target, proposals can make large efficient moves and cross separated modes; $g=\pi$ gives acceptance one and independent samples. On the other hand, a good global approximation is required. A proposal with tails too light can leave the chain trapped at states having very large importance weight $\pi /g$, since most outgoing proposals then have tiny acceptance probabilities.

For perspective, if the normalized target satisfies $\pi \le Mg$, the accepted transition density is at least $\pi ({\theta }^{\prime })/M$, because both entries in $\min \{g({\theta }^{\prime }),\pi ({\theta }^{\prime })/w(\theta )\}$ have that lower bound, where $w=\pi /g\le M$. This yields a uniform refresh component. Thus a well-designed independence proposal can mix rapidly, whereas a poor global proposal can be much worse than a locally tuned random walk. The tradeoff concerns effective samples per computation, not acceptance rate alone.

Let $v={\sigma }^2>0$, $d_k=k+1$, and $S_k(a)=\Vert y-X_{k}a{\Vert }^2$. Under the standard interpretation that the two stated priors are independent, multiply the normal likelihood by the normal coefficient prior and the inverse-gamma distribution density. The independent normal and inverse-gamma regression priors giveFor a fixed model the determinant and Gaussian normalizing constants are independent of $a_k,v$ and may be omitted here; ${\Sigma }_k$ is assumed positive-definite. The coefficient prior is not scaled by $v$, so replacing its precision by ${\Sigma }_k^{-1}/v$ would describe a different prior.

For an additional check, completing the coefficient square and reading the variance kernel give the full conditional distributionsandMarginal priors alone would not determine their joint prior without the independence assumption. The front-page inverse-gamma mean has a typographical error: this density has mean $\beta /(\alpha -1)$ when $\alpha >1$, not a denominator $\alpha +1$. The density itself, used in the posterior calculation, fixes the correct parameterization.

The cross-model target must include a model-order prior ${\rho }_k$. Write its joint, unnormalized density asHere the prior densities are normalized within each model. In particular their $k$-dependent determinants and powers of $2\pi$ cannot be dropped merely because the fixed-model posterior in part (i) was given up to proportionality.

For explicit birth and death moves for Bayesian variable selection, choose a birth with probability $b_k$, generate $u$ from a density $g_k(u\mid a,v)$, append $u$ as the new covariate coefficient, and leave $v$ unchanged. The map from $(a,v,u)$ to $(a^{\prime },v^{\prime })=((a,u),v)$ is invertible under the reverse deletion and has absolute Jacobian one. It matches dimensions: the smaller-model coefficients and variance have dimension $k+2$, and the one auxiliary draw supplies the larger-model dimension $k+3$. If the reverse death is selected with probability $d_{k+1}$, accept withFor a death, delete the last coefficient $u$, retain the others and the variance, and useTo verify detailed balance, multiply the birth acceptance by $t_k{b}_k{g}_k$ and the reverse death acceptance by $t_{k+1}{d}_{k+1}$. Both equal the minimum of these two quantities in the matched coordinates. This proves reversibility without merely naming a dimension-changing algorithm. A more general bijection introduces its absolute Jacobian, and random choices among several covariates introduce the corresponding forward and reverse selection probabilities. Add within-model updates so that coefficients and variance can explore each model. These are reversible-jump Markov chain Monte Carlo updates; the unspecified model prior and proposal density must be supplied to implement them numerically.

First distinguish a complete-data calculation, conditional on the missing counts, from optimization of the observed-data likelihood. Introduce death counts $D_j={\sum }_{i=1}^j{n}_{ij}$, survivor counts $S_j=A_{1j}$, observed recapture counts $C_j={\sum }_{i=1}^{j-1}{m}_{ij}$ and unobserved-alive counts $U_j=A_{2j}$. Terms involving these probabilities in the complete-data log likelihood areThe survival score is $S_j/{\varphi }_j-D_j/(1-{\varphi }_j)$ and the detection score is $C_j/P_j-U_j/(1-P_j)$. Their derivatives are nonpositive, and strictly negative whenever their respective total counts are positive. The interior roots, with boundary interpretations if a success or failure count vanishes, are thereforeSubstituting the survivor and unobserved-alive totals givesThese are the requested complete-data maximizing ratios. Here $n_{iJ}$ denotes the terminal missing count discussed in part (b); pre-release cells and nonexistent release cohorts contribute zero. If a denominator is zero, the likelihood is constant in that parameter rather than defining a $0/0$ estimate. Also $P_1$ does not occur in this likelihood and has no identified maximum. With unobserved $n$, these ratios alone do not yet supply an observed-data estimate; the EM algorithm below supplies the missing-count expectations.

The consistent missing-data interpretation is that each never-recaptured individual either dies in an interval $j<J$ or survives unobserved through the final occasion. Define the known never-recaptured total and its allocation byThe terminal cell $n_{iJ}$ is the source's separately named $n_i$. The printed identity equating that cell to the total observed recaptures is inconsistent with the supplied multinomial law. In particular, with no observed recaptures its printed value would be zero even though the terminal-cell probability is positive at interior parameters. The displayed accounting identity repairs that notation and makes all sums through $J$ well-defined.

For an individual released at $i$, define the unnormalized probabilitiesThe empty product is one. Each intermediate factor means surviving an interval and avoiding its subsequent detection; the death weight additionally includes failure to survive interval $j$. Let $K_i={\sum }_{j=i}^J{w}_{ij}$. This is the probability of never being recaptured. It is the source's normalizing constant $C_i$, renamed to avoid confusion with the detection-success counts above. Conditional on no recapture, division by $K_i$ gives precisely the supplied multinomial distribution cell probabilities.

Initialize an interior parameter vector $({\varphi }^{(0)},P^{(0)})$. At iteration $r$, the E step of missing-count EM for capture-recapture computes, for every cohort and missing cell,If $N_i=0$, set every missing count to zero without dividing, even if $K_i=0$. If $N_i>0$ and $K_i=0$, that parameter vector gives zero probability to the observations and cannot be used in the E step. These expected counts sum to $N_i$, preserving each cohort's observed accounting. The E step is conditional expectation, not replacement by the modal missing counts or a single random imputation.

The expected complete-data log likelihood is linear in the missing counts, apart from count-only terms irrelevant to the new parameter values. In the M step substitute these expectations into the success/failure ratios from part (a):andIf an exposure denominator vanishes, retain any admissible value of the unidentifiable parameter. Repeat until the observed likelihood and parameter updates stabilize. The conditional normalizers in the E step are evaluated at the old parameters; maximizing a newly normalized conditional missing-data distribution instead of the expected complete-data log likelihood would be the wrong M step.

One can monitor the observed likelihood directly. The probability of a first recapture at $s>i$ isand successive death, detection and continuation alternatives partition the possibilities, giving $K_i+{\sum }_{s=i+1}^J{p}_{is}=1$. Thus, up to count-only factors,The usual conditional Jensen argument for the EM algorithm makes this likelihood nondecreasing at each exact M step. It does not guarantee a unique global maximum. In this observation scheme the final survival and detection probabilities enter through the product ${\varphi }_{J-1}{P}_J$ only, so these two probabilities cannot generally be estimated separately without a further constraint. The algorithm remains valid, but may converge to different points on the same likelihood ridge.

With complete counts fixed, the likelihood factor for $P_j$ is $P_j^{C_j}(1-P_j{)}^{U_j}$. As a function of $P_j\in (0,1)$, normalize it to the density of a beta distribution with parameters $C_j+1,U_j+1$. Normalization does not change its maximizing argument. The standard beta-mode formula, when both counts are positive, givesThis uses the standard distribution result without differentiating the likelihood, and introducing the normalized beta kernel does not require adopting a Bayesian prior. If only $C_j$ is zero the maximum is at zero; if only $U_j$ is zero it is at one; if both vanish every probability maximizes the constant factor. For the EM M step the same beta-mode argument applies to the expected, possibly noninteger counts, since positive real beta parameters are allowed.