1\section{Metal execution and precision strategy}2\label{sec:metal}34\paragraph{The constraint.} Apple GPUs and the Metal Shading Language have no 64-bit floating-point type; a kernel declaring a \texttt{double} does not compile (``\texttt{'double' is not supported in Metal}''). Econometric accumulators---sums of $4\times10^{5}$ squared returns of magnitude $10^{-6}$, Gram matrices $\bX'\bX$ with regressors of very different scales, log-likelihoods summed over thousands of observations---lose several digits in naive FP32, whose unit roundoff is $u = 2^{-24} \approx 6\times10^{-8}$. The strategy has four parts.56\paragraph{1. Double-float accumulation.} Every accumulator on the GPU is a pair $(h, \ell)$ of \texttt{float}s representing $h + \ell$ with $\abs{\ell} \le \tfrac12\ulp(h)$. Additions use the error-free transformation TwoSum \citep{Dekker1971,Neumaier1974}: for floats $a, b$,7\begin{equation}8s = \fl(a + b), \qquad e = \bigl(a - (s - b')\bigr) + \bigl(b - b'\bigr),\quad b' = \fl(s - a), \qquad a + b = s + e \ \text{exactly},9\label{eq:twosum}10\end{equation}11followed by a renormalisation of the pair. The pair carries about 48 significant bits, and the analysis of \citet{Higham1993} bounds the relative error of a compensated sum of $n$ terms by $2u + O(nu^{2})$ instead of $(n-1)u$ for recursive summation. The Metal compiler must not reassociate: the kernels are compiled with \texttt{-fno-fast-math} offline and \texttt{fastMathEnabled = NO} at runtime. The inputs are still rounded to FP32 on the way in; for a sum of squares of $M$ independent returns with rounding errors $\delta_j$, $\abs{\delta_j} \le u$,12\begin{equation}13\frac{\fl(\RV) - \RV}{\RV} \;\approx\; \frac{2\sum_j \delta_j\, r_j^{2}}{\sum_j r_j^{2}}, \qquad14\sqrt{\E\Bigl[\Bigl(\tfrac{2\sum_j \delta_j r_j^{2}}{\sum_j r_j^{2}}\Bigr)^{2}\Bigr]} \;\lesssim\; 2u\,\sqrt{\frac{\sum_j r_j^{4}}{\bigl(\sum_j r_j^{2}\bigr)^{2}}} \;\approx\; 2u\sqrt{\frac{3}{M}}15\label{eq:rounding}16\end{equation}17for Gaussian returns, i.e. about $10^{-8}$ for $M = 390$ and decreasing in $M$. Figure~\ref{fig:precision} shows that the maximum error observed across the segments of each E1 case follows this $M^{-1/2}$ law and lies within a factor of about two of the expected maximum of 926 Gaussian deviations ($\approx 3.2$ root-mean-square units), the excess at short segments reflecting the heavy tails of real returns.1819\begin{figure}[t]20\centering\includegraphics[width=.62\textwidth]{figures/fig_precision.pdf}21\caption{Precision of the Metal realized variance. Maximum relative difference between the Metal (FP32 double-float) and the CPU FP64 realized variance across the segments of each E1 case, versus the number of returns per segment, with the root-mean-square prediction of \eqref{eq:rounding} and its expected maximum over 926 segments.}22\label{fig:precision}23\end{figure}2425\paragraph{2. Chunked partial sums combined in FP64.} A segment longer than the chunk length ($4\,096$ elements by default, a \texttt{PrecisionPolicy} parameter) is split into chunks; one threadgroup of 256 threads reduces each chunk with a fixed-order tree reduction of double-float values and writes a \texttt{float2}; the host adds the chunks in FP64 in a fixed order and adds the chunk-boundary cross terms ($\abs{r_j}\abs{r_{j-1}}$ across a boundary) in FP64. The combination order is deterministic, so results are bitwise reproducible across runs for a given chunk length and threadgroup width (Section~\ref{sec:correct}).2627\paragraph{3. Mixed precision by design.} Quantities that are cheap and precision-critical stay in FP64 on the CPU: the finite-sample constants, the initial variance of a GARCH recursion, the Student-$t$ constant $c(\nu)$ and $\E\abs{z}$ in \eqref{eq:garch-ll}--\eqref{eq:egarch} (Metal has no \texttt{lgamma}), the regression solve from the GPU-accumulated sufficient statistics, and all inference. The batched GARCH kernel scales each series to unit variance on the host, $r_t' = r_t/s$, $\mu' = \mu/s$, $\omega' = \omega/s^{2}$ (EGARCH: $\omega' = \omega + (1-\beta)\log s^{-2}$), and the host un-scales the log-likelihood by $-T\log s$, which conditions the FP32 recursion for returns of magnitude $10^{-2}$.2829\paragraph{4. Kernel fusion.} The realized kernel computes the five accumulators of \eqref{eq:rv-bpv}--\eqref{eq:rq-j} in a single pass over the data; the four-kernel variant exists only for the ablation. Regression statistics accumulate the upper triangle of $\bX'\bX$, $\bX'\by$ and $\by'\by$ ($k \le 8$) in registers. The GARCH kernel runs one recursion per thread for one (series, parameter vector) pair---the parallelism comes from the number of parameter vectors, not from the recursion.3031\paragraph{Host side.} The backend owns one device, one command queue, one library and one pipeline cache per process; buffers come from a pool keyed by power-of-two size class in shared storage (unified memory). The shader library is loaded from a prebuilt \texttt{.metallib} when present and compiled from the embedded MSL source otherwise; both paths are measured. GPU time is read from the command buffer's \texttt{GPUStartTime}/\texttt{GPUEndTime}; wall-clock includes FP64$\to$FP32 conversion into the shared buffer, encoding, commit and completion. Metal does not expose the GPU core count; the manifest records it from \texttt{system\_profiler}.3233\paragraph{What the GPU does not do.} Optimisation drivers (BFGS line searches) run on the CPU; the GPU evaluates batches of likelihoods. Rolling sums use a direct $O(w)$-per-output kernel, which the experiments show to be slower than the CPU's $O(1)$ running sum---a limitation we report rather than hide. Nothing in the Metal path is a new estimator.34