[proofplan] Writing $Y_t = \sigma_t^2$ and $A_t = \alpha Z_{t-1}^2 + \beta$, the volatility recursion becomes the affine random recurrence $Y_t = \omega + A_t Y_{t-1}$. Iterating it backward exhibits the backward products $B_{t,k} = \prod_{j=1}^k(\alpha Z_{t-j}^2 + \beta)$, and the heart of the proof is that the Nelson condition $\mathbb E[\log(\alpha Z_0^2 + \beta)] < 0$ forces these products to decay geometrically almost surely, via the Strong Law of Large Numbers applied (after a truncation that handles a possibly infinite log-moment) to the i.i.d. summands $\log(\alpha Z_{t-j}^2 + \beta)$. Geometric decay makes the series $\omega \sum_k B_{t,k}$ converge almost surely; a rearrangement of its nonnegative terms shows it solves the recursion, and it is manifestly a fixed measurable function of the past innovations, hence causal and strictly stationary. Uniqueness follows because the difference of any two strictly stationary causal solutions equals $B_{t,n}$ times a remote difference, which vanishes as $n \to \infty$ — almost surely for the explicit solution and in probability for the competitor, using only that its marginal law is fixed. [/proofplan]

[step:Reduce the volatility recursion to an affine random recurrence] Define the random coefficients and the volatility sequence \begin{align*} A_t &:= \alpha Z_{t-1}^2 + \beta, & Y_t &:= \sigma_t^2 \qquad (t \in \mathbb Z), \end{align*} so that each $A_t \ge 0$ is $\sigma(Z_{t-1})$-measurable and the family $(A_t)_{t \in \mathbb Z}$ is i.i.d. with the law of $A_0 = \alpha Z_0^2 + \beta$. Substituting $X_{t-1} = \sigma_{t-1} Z_{t-1}$ into the defining recursion gives \begin{align*} \sigma_t^2 = \omega + \alpha (\sigma_{t-1} Z_{t-1})^2 + \beta \sigma_{t-1}^2 = \omega + (\alpha Z_{t-1}^2 + \beta)\sigma_{t-1}^2, \end{align*} that is, \begin{align*} Y_t = \omega + A_t Y_{t-1} \qquad (t \in \mathbb Z). \tag{R} \end{align*} Iterating (R) backward $n$ times yields, for every $n \ge 1$, \begin{align*} Y_t = \omega \sum_{k=0}^{n-1} \Big(\prod_{i=0}^{k-1} A_{t-i}\Big) + \Big(\prod_{i=0}^{n-1} A_{t-i}\Big) Y_{t-n}. \tag{R$_n$} \end{align*} With the substitution $i = j-1$ we have $\prod_{i=0}^{k-1} A_{t-i} = \prod_{j=1}^{k} A_{t-j+1} = \prod_{j=1}^{k}(\alpha Z_{t-j}^2 + \beta)$, since $A_{t-j+1} = \alpha Z_{(t-j+1)-1}^2 + \beta = \alpha Z_{t-j}^2 + \beta$. We therefore set \begin{align*} B_{t,k} := \prod_{j=1}^{k}\big(\alpha Z_{t-j}^2 + \beta\big), \qquad B_{t,0} := 1, \end{align*} so that $B_{t,k}$ is a nonnegative $\mathcal F_{t-1}$-measurable random variable, finite $\mathbb P$-a.s. (a finite product of a.s. finite factors), and (R$_n$) reads $Y_t = \omega \sum_{k=0}^{n-1} B_{t,k} + B_{t,n} Y_{t-n}$. [/step]

[step:Prove almost sure geometric decay of the backward products]Fix $t \in \mathbb Z$. For $j \ge 1$ set $\xi_j := \log\big(\alpha Z_{t-j}^2 + \beta\big) \in [-\infty, \infty)$, so that the $(\xi_j)_{j \ge 1}$ are i.i.d. with the law of $\xi_0 := \log(\alpha Z_0^2 + \beta)$, and $\log B_{t,k} = \sum_{j=1}^{k} \xi_j$ (with the convention $\log 0 = -\infty$). We claim there is a constant $c < 0$ such that \begin{align*} \limsup_{k \to \infty} \frac{1}{k}\log B_{t,k} = \limsup_{k \to \infty} \frac{1}{k}\sum_{j=1}^{k}\xi_j \le c < 0 \qquad \mathbb P\text{-a.s.} \tag{D} \end{align*} Choose any $c \in (\gamma, 0)$, which is possible because $\gamma < 0$. For $M > 0$ define the truncations \begin{align*} \xi_j^{(M)} := \max(\xi_j, -M), \end{align*} which are i.i.d. and integrable: $0 \le (\xi_j^{(M)})^+ = \xi_j^+ \le \alpha Z_{t-j}^2 + \beta$ has expectation $\le \alpha + \beta < \infty$, while $(\xi_j^{(M)})^- \le M$. Writing $\xi_0^{(M)} = \xi_0^+ - (\xi_0^- \wedge M)$ and noting $\xi_0^- \wedge M \uparrow \xi_0^-$ as $M \uparrow \infty$, the [Monotone Convergence Theorem](/theorems/509) gives $\mathbb E[\xi_0^- \wedge M] \uparrow \mathbb E[\xi_0^-]$, hence \begin{align*} \gamma_M := \mathbb E\big[\xi_0^{(M)}\big] = \mathbb E[\xi_0^+] - \mathbb E[\xi_0^- \wedge M] \;\downarrow\; \mathbb E[\xi_0^+] - \mathbb E[\xi_0^-] = \gamma. \end{align*} Since $\gamma_M \downarrow \gamma < c$, we may fix $M$ with $\gamma_M < c$. The $\xi_j^{(M)}$ are i.i.d. and integrable, so the [Strong Law of Large Numbers](/theorems/520) yields $\frac1k \sum_{j=1}^k \xi_j^{(M)} \to \gamma_M$ $\mathbb P$-a.s. As $\xi_j \le \xi_j^{(M)}$ pointwise, \begin{align*} \limsup_{k \to \infty} \frac{1}{k}\sum_{j=1}^{k}\xi_j \le \lim_{k \to \infty} \frac{1}{k}\sum_{j=1}^{k}\xi_j^{(M)} = \gamma_M < c \qquad \mathbb P\text{-a.s.}, \end{align*} which is (D). In particular $\frac1k \log B_{t,k} \to$ a value $\le c < 0$, so $\log B_{t,k} \to -\infty$ and \begin{align*} B_{t,n} \to 0 \qquad \mathbb P\text{-a.s. as } n \to \infty. \tag{D$'$} \end{align*}[/step]

[guided]We must turn the qualitative condition $\gamma < 0$ into a quantitative, pathwise decay rate for the products $B_{t,k} = \prod_{j=1}^k(\alpha Z_{t-j}^2 + \beta)$. The natural device is to take logarithms, converting the product into a sum of i.i.d. terms $\xi_j = \log(\alpha Z_{t-j}^2 + \beta)$ to which a law of large numbers applies. The $(\xi_j)_{j\ge1}$ are i.i.d. because they are measurable functions of the distinct innovations $Z_{t-1}, Z_{t-2}, \dots$, and they share the law of $\xi_0 = \log(\alpha Z_0^2 + \beta)$. Why can we not apply the Strong Law directly to $(\xi_j)$? The classical statement requires integrability, i.e. $\mathbb E|\xi_0| < \infty$. Here the **positive** part is controlled — using $\log^+ x \le x$ for $x \ge 0$ we get $\xi_0^+ = \log^+(\alpha Z_0^2 + \beta) \le \alpha Z_0^2 + \beta$, whose expectation is $\alpha\,\mathbb E[Z_0^2] + \beta = \alpha + \beta < \infty$ thanks to the normalisation $\mathbb E[Z_0^2] = 1$. But the **negative** part may be infinite: if $\beta = 0$ and $\mathbb P(Z_0 = 0) > 0$, then $\xi_0 = -\infty$ on a set of positive probability, so $\mathbb E[\xi_0^-] = \infty$ and $\gamma = -\infty$. Both regimes ($\gamma$ finite and $\gamma = -\infty$) must be handled. The remedy is truncation from below. Replace $\xi_j$ by $\xi_j^{(M)} = \max(\xi_j, -M)$. These are bounded below by $-M$ and retain the integrable positive part, hence are integrable, so the [Strong Law of Large Numbers](/theorems/520) applies and gives $\frac1k\sum_{j=1}^k \xi_j^{(M)} \to \gamma_M := \mathbb E[\xi_0^{(M)}]$ almost surely. How is $\gamma_M$ related to $\gamma$? Decomposing $\xi_0^{(M)} = \xi_0^+ - (\xi_0^- \wedge M)$, the truncated negative part $\xi_0^- \wedge M$ increases to $\xi_0^-$ as $M \to \infty$, so by the [Monotone Convergence Theorem](/theorems/509) its expectation increases to $\mathbb E[\xi_0^-]$, and therefore $\gamma_M$ decreases to $\gamma$. Since $\gamma < 0$, pick any target level $c$ strictly between $\gamma$ and $0$; because $\gamma_M \downarrow \gamma < c$, some finite truncation level $M$ already satisfies $\gamma_M < c$. Now we exploit the one-sided inequality $\xi_j \le \xi_j^{(M)}$ (truncation only raises values). Summing and dividing by $k$, \begin{align*} \limsup_{k\to\infty} \frac1k \sum_{j=1}^k \xi_j \le \lim_{k\to\infty}\frac1k\sum_{j=1}^k \xi_j^{(M)} = \gamma_M < c < 0 \qquad \mathbb P\text{-a.s.} \end{align*} This is exactly (D). Geometrically, the average log-factor is eventually below the strictly negative level $c$, so $\log B_{t,k}$ drifts to $-\infty$ and the products collapse to $0$ — statement (D$'$). This collapse is the engine of both existence (a convergent series) and uniqueness (a vanishing remote term).[/guided]

[step:Define the explicit solution and verify convergence, positivity, and causality] For each $t \in \mathbb Z$ define \begin{align*} \sigma_t^2 := \omega \sum_{k=0}^{\infty} B_{t,k} = \omega \sum_{k=0}^{\infty} \prod_{j=1}^{k}\big(\alpha Z_{t-j}^2 + \beta\big). \end{align*} By (D), $\mathbb P$-a.s. there is a (random) index $K$ with $B_{t,k} \le e^{ck}$ for all $k \ge K$. Splitting the series, \begin{align*} \sum_{k=0}^{\infty} B_{t,k} \le \sum_{k=0}^{K-1} B_{t,k} + \sum_{k=K}^{\infty} e^{ck} < \infty \qquad \mathbb P\text{-a.s.}, \end{align*} since the head is a finite sum of a.s. finite terms and the tail is a convergent geometric series ($c < 0$). Thus $\sigma_t^2$ is finite $\mathbb P$-a.s.; it is nonnegative, and in fact $\sigma_t^2 \ge \omega B_{t,0} = \omega > 0$. As a series of $\mathcal F_{t-1}$-measurable functions, $\sigma_t^2$ is $\mathcal F_{t-1}$-measurable, hence the solution is causal. It remains to verify the recursion (R). All terms below are nonnegative, so regrouping the series is justified ([Fubini–Tonelli Theorem](/theorems/513) for the counting measure on $\mathbb N_0$). Using $A_t = \alpha Z_{t-1}^2 + \beta$ and the identity $A_t\,B_{t-1,k} = A_t \prod_{j=1}^{k}(\alpha Z_{t-1-j}^2 + \beta) = \prod_{i=1}^{k+1}(\alpha Z_{t-i}^2 + \beta) = B_{t,k+1}$, \begin{align*} \omega + (\alpha Z_{t-1}^2 + \beta)\,\sigma_{t-1}^2 &= \omega + A_t\,\omega\sum_{k=0}^{\infty} B_{t-1,k} = \omega + \omega\sum_{k=0}^{\infty} B_{t,k+1} \\ &= \omega\sum_{m=0}^{\infty} B_{t,m} = \sigma_t^2, \end{align*} where the last sum reindexes $m = k+1$ and absorbs $\omega = \omega B_{t,0}$. Hence $(\sigma_t^2)_{t \in \mathbb Z}$ satisfies (R) for every $t$. [/step]

[step:Show the explicit solution is strictly stationary]Define the measurable map \begin{align*} g : \mathbb R^{\mathbb N} &\to [0, \infty] \\ (z_1, z_2, \dots) &\mapsto \omega \sum_{k=0}^{\infty} \prod_{j=1}^{k}\big(\alpha z_j^2 + \beta\big), \end{align*} which is measurable as a pointwise limit of the measurable partial sums (each partial sum is a polynomial in finitely many coordinates). By construction $\sigma_t^2 = g(Z_{t-1}, Z_{t-2}, \dots)$ for every $t$, and $X_t = \sigma_t Z_t = \sqrt{g(Z_{t-1}, Z_{t-2}, \dots)}\,Z_t$, so each $\sigma_t^2$ and each $X_t$ is a fixed measurable function of the innovations indexed at or before $t$, the same function for every $t$. Fix a finite collection $t_1 < \dots < t_m$ and $\tau \in \mathbb Z$. Since $\sigma_{t_i}^2$ depends only on $Z_s$ with $s \le t_i - 1 \le t_m - 1$, there is a single measurable map $F : \mathbb R^{\mathbb N_0} \to \mathbb R^m$, independent of $\tau$, with \begin{align*} (\sigma_{t_1}^2, \dots, \sigma_{t_m}^2) &= F\big((Z_{t_m - i})_{i \ge 0}\big), & (\sigma_{t_1 + \tau}^2, \dots, \sigma_{t_m + \tau}^2) &= F\big((Z_{t_m + \tau - i})_{i \ge 0}\big). \end{align*} The two backward sequences $(Z_{t_m - i})_{i \ge 0}$ and $(Z_{t_m + \tau - i})_{i \ge 0}$ are both i.i.d. sequences with the common marginal law of $Z_0$, so they have identical laws on $\mathbb R^{\mathbb N_0}$. Applying the same measurable $F$ to equal-in-law inputs produces equal-in-law outputs: \begin{align*} (\sigma_{t_1 + \tau}^2, \dots, \sigma_{t_m + \tau}^2) \overset{d}{=} (\sigma_{t_1}^2, \dots, \sigma_{t_m}^2). \end{align*} As this holds for every finite collection and every shift $\tau$, the process $(\sigma_t^2)_{t \in \mathbb Z}$ is strictly stationary. Together with Step 3 it is a strictly stationary causal solution, establishing existence.[/step]

[guided]Strict stationarity asks that all finite-dimensional distributions are invariant under time shifts. The clean way to see this is to recognise $(\sigma_t^2)$ as a **factor of an i.i.d. sequence**: each coordinate is the *same* deterministic measurable recipe $g$ applied to the past innovations. Two features make the argument work. First, $g$ does not depend on $t$ — this is forced by the autonomous form of the recursion and is exactly why we wrote the solution as a function of the relative lags $Z_{t-1}, Z_{t-2}, \dots$ rather than of absolute time. Second, the input process $(Z_t)_{t\in\mathbb Z}$ is i.i.d., hence its law is invariant under the shift $\tau$. Concretely, bundle the finitely many times $t_1 < \dots < t_m$. Each $\sigma_{t_i}^2$ reads off innovations no later than $t_m - 1$, so the whole vector is one measurable function $F$ of the backward stream $(Z_{t_m - i})_{i \ge 0}$. Shifting all times by $\tau$ replaces this stream by $(Z_{t_m + \tau - i})_{i \ge 0}$ but keeps the *same* $F$. Why do the two streams have the same law? Each is an i.i.d. sequence drawn from the law of $Z_0$; the absolute starting index is irrelevant to the joint distribution of an i.i.d. family. Pushing equal laws through the common measurable map $F$ gives equal laws of the output vectors — precisely shift-invariance of the finite-dimensional distributions. What would fail without the i.i.d. (or at least stationarity) assumption on $(Z_t)$? The two input streams could have different laws, and the conclusion would collapse.[/guided]

[step:Prove uniqueness up to indistinguishability]Let $(\hat\sigma_t^2)_{t \in \mathbb Z}$ be any strictly stationary causal solution, and let $(\sigma_t^2)$ be the explicit solution of Steps 3–4. Both satisfy (R$_n$): for every $n \ge 1$, \begin{align*} \sigma_t^2 &= \omega\sum_{k=0}^{n-1} B_{t,k} + B_{t,n}\,\sigma_{t-n}^2, & \hat\sigma_t^2 &= \omega\sum_{k=0}^{n-1} B_{t,k} + B_{t,n}\,\hat\sigma_{t-n}^2. \end{align*} Subtracting cancels the common partial sum: \begin{align*} \sigma_t^2 - \hat\sigma_t^2 = B_{t,n}\big(\sigma_{t-n}^2 - \hat\sigma_{t-n}^2\big) \qquad (n \ge 1). \tag{$\ast$} \end{align*} The left-hand side does not depend on $n$; we show the right-hand side tends to $0$ in probability. For the explicit solution, the identity $B_{t,n} B_{t-n,k} = B_{t,n+k}$ (the same reindexing as in Step 3) gives \begin{align*} B_{t,n}\,\sigma_{t-n}^2 = \omega \sum_{k=0}^{\infty} B_{t,n} B_{t-n,k} = \omega\sum_{k=0}^{\infty} B_{t,n+k} = \omega\sum_{m=n}^{\infty} B_{t,m} \xrightarrow{n\to\infty} 0 \quad \mathbb P\text{-a.s.}, \end{align*} the tail of the $\mathbb P$-a.s. convergent series from Step 3. For the competitor, fix $\varepsilon > 0$ and $M > 0$. On the event $\{\hat\sigma_{t-n}^2 \le M\}$ we have $B_{t,n}\hat\sigma_{t-n}^2 > \varepsilon \implies B_{t,n} > \varepsilon/M$, so \begin{align*} \mathbb P\big(B_{t,n}\hat\sigma_{t-n}^2 > \varepsilon\big) \le \mathbb P\big(B_{t,n} > \varepsilon/M\big) + \mathbb P\big(\hat\sigma_{t-n}^2 > M\big). \end{align*} By strict stationarity $\hat\sigma_{t-n}^2 \overset{d}{=} \hat\sigma_0^2$, so the second term equals $\mathbb P(\hat\sigma_0^2 > M)$. By (D$'$), $B_{t,n} \to 0$ $\mathbb P$-a.s. and hence in probability, so the first term $\to 0$ as $n \to \infty$ for fixed $M$. Therefore \begin{align*} \limsup_{n \to \infty} \mathbb P\big(B_{t,n}\hat\sigma_{t-n}^2 > \varepsilon\big) \le \mathbb P\big(\hat\sigma_0^2 > M\big), \end{align*} and letting $M \to \infty$ — using that $\hat\sigma_0^2$ is $\mathbb P$-a.s. finite — gives $\mathbb P(\hat\sigma_0^2 > M) \to 0$. Hence $B_{t,n}\hat\sigma_{t-n}^2 \to 0$ in probability. Consequently the right-hand side of $(\ast)$, namely $B_{t,n}\sigma_{t-n}^2 - B_{t,n}\hat\sigma_{t-n}^2$, converges to $0$ in probability. Its left-hand side $\sigma_t^2 - \hat\sigma_t^2$ is constant in $n$, so it must equal its in-probability limit: $\sigma_t^2 = \hat\sigma_t^2$ $\mathbb P$-a.s., for each fixed $t$. Since $\mathbb Z$ is countable, the exceptional null sets unite to a single null set, so $\mathbb P(\sigma_t^2 = \hat\sigma_t^2 \text{ for all } t) = 1$; the two solutions are indistinguishable.[/step]

[guided]Why should two strictly stationary causal solutions coincide? Both are pinned to the *same* finite expansion in (R$_n$): iterating the recursion $n$ steps backward produces the identical partial sum $\omega\sum_{k=0}^{n-1}B_{t,k}$ for either solution, plus a remote remainder $B_{t,n}$ times a value $n$ steps in the past. Subtracting kills the shared partial sum and leaves identity $(\ast)$: the present discrepancy equals the backward product $B_{t,n}$ times a remote discrepancy. The whole game is to send the remote remainders to zero. The crucial observation is that the left side of $(\ast)$ is the *same* number for every $n$ — it is a fixed random variable — so if the right side has any limit (in any reasonable mode), the left side equals that limit. For the explicit solution this is transparent: $B_{t,n}\sigma_{t-n}^2$ telescopes via $B_{t,n}B_{t-n,k} = B_{t,n+k}$ into the tail $\omega\sum_{m \ge n} B_{t,m}$ of an almost surely convergent series, which vanishes almost surely. The competitor $\hat\sigma_{t-n}^2$ is the delicate term: we cannot bound it pathwise, because we know nothing about it beyond the three defining properties. But we know two soft facts. First, $B_{t,n}$ — built only from $Z_{t-1}, \dots, Z_{t-n}$ — tends to $0$ almost surely by Step 2. Second, strict stationarity pins the *marginal law* of $\hat\sigma_{t-n}^2$ to that of $\hat\sigma_0^2$, a fixed almost surely finite random variable, so it cannot escape to infinity in distribution. We convert this into convergence in probability by a truncation at level $M$: either $\hat\sigma_{t-n}^2$ exceeds $M$ (an event of fixed small probability $\mathbb P(\hat\sigma_0^2 > M)$, by stationarity), or it is at most $M$, in which case $B_{t,n}\hat\sigma_{t-n}^2 > \varepsilon$ requires $B_{t,n} > \varepsilon/M$, which becomes improbable as $n \to \infty$ because $B_{t,n} \to 0$. Taking $n \to \infty$ and then $M \to \infty$ drives the probability to $0$. Notice we used only that $\hat\sigma_{t-n}^2$ has a fixed marginal law and is finite — exactly what strict stationarity and the solution property supply; no independence between $B_{t,n}$ and $\hat\sigma_{t-n}^2$ is needed. With both remote terms vanishing (one a.s., one in probability), the right side of $(\ast)$ tends to $0$ in probability, forcing the constant left side $\sigma_t^2 - \hat\sigma_t^2$ to be $0$ almost surely for each $t$. Finiteness of the index set $\mathbb Z$ up to countability then upgrades "almost surely for each $t$" to "almost surely for all $t$ simultaneously," which is indistinguishability. Combined with the existence from Steps 3–4, this proves the theorem.[/guided]