Prove that F-statistic follows F-distribution

In light of this question : Proof that the coefficients in an OLS model follow a t-distribution with (n-k) degrees of freedom

I would love to understand why

$$F = \frac{(\text{TSS}-\text{RSS})/(p-1)}{\text{RSS}/(n-p)},$$

where $p$ is the number of model parameters and $n$ the number of observations and $TSS$ the total variance, $RSS$ the residual variance, follows an $F_{p-1,n-p}$ distribution.

I must admit I have not even attempted to prove it as I wouldn't know where to start.

• Christoph Hanck and Francis has given a very good answer already. If you still have difficulties in understanding the proof of f test for linear regression, try to checkout teamdable.github.io/techblog/… . I wrote the blog post about the proof of the ftest for linear regression. It is written in Korean but it may not be a problem because almost all of it is math formula. I hope it would help if you still have difficulties in understanding the proof of f test for linear regression. Commented Aug 5, 2019 at 1:34
• While this link may answer the question, it is better to include the essential parts of the answer here and provide the link for reference. Link-only answers can become invalid if the linked page changes. - From Review
– mkt
Commented Aug 5, 2019 at 5:36
• Ultimately, there are only three fundamental things to know here. The first is that sums of squares of zero-mean Normal variables are multiples of chi-squared distributions. This is often taken as the definition of a chi-square distribution. The second is that TSS-RSS and RSS are independent; this is a matter of linear algebra, which expresses them as functions of uncorrelated Normal variables. The third is that a ratio of independent chi-squared variables has (up to a constant multiple) an $F$ distribution: indeed, you can take this as a definition of $F.$
– whuber
Commented Feb 4, 2022 at 16:11
• It must be stressed that $F$ follows $F_{p - 1, n - p}$ distribution only when the null hypothesis $\beta_1 = \cdots = \beta_{p - 1} = 0$ holds. In other words, without this condition, it does not have $F$ distribution generally. Commented Mar 15 at 10:54

Let us show the result for the general case of which your formula for the test statistic is a special case. In general, we need to verify that the statistic can be, according to the characterization of the $$F$$ distribution, be written as the ratio of independent $$\chi^2$$ r.v.s divided by their degrees of freedom.

Let $$H_{0}:R^\prime\beta=r$$ with $$R$$ and $$r$$ known, nonrandom and $$R:k\times q$$ has full column rank $$q$$. This represents $$q$$ linear restrictions for (unlike in OPs notation) $$k$$ regressors including the constant term. So, in @user1627466's example, $$p-1$$ corresponds to the $$q=k-1$$ restrictions of setting all slope coefficients to zero.

In view of $$Var\bigl(\hat{\beta}_{\text{ols}}\bigr)=\sigma^2(X'X)^{-1}$$, we have $$\begin{eqnarray*} R^\prime(\hat{\beta}_{\text{ols}}-\beta)\sim N\left(0,\sigma^{2}R^\prime(X^\prime X)^{-1} R\right), \end{eqnarray*}$$ so that (with $$B^{-1/2}=\{R^\prime(X^\prime X)^{-1} R\}^{-1/2}$$ being a "matrix square root" of $$B^{-1}=\{R^\prime(X^\prime X)^{-1} R\}^{-1}$$, via, e.g., a Cholesky decomposition) $$\begin{eqnarray*} n:=\frac{B^{-1/2}}{\sigma}R^\prime(\hat{\beta}_{\text{ols}}-\beta)\sim N(0,I_{q}), \end{eqnarray*}$$ as $$\begin{eqnarray*} Var(n)&=&\frac{B^{-1/2}}{\sigma}R^\prime Var\bigl(\hat{\beta}_{\text{ols}}\bigr)R\frac{B^{-1/2}}{\sigma}\\ &=&\frac{B^{-1/2}}{\sigma}\sigma^2B\frac{B^{-1/2}}{\sigma}=I \end{eqnarray*}$$ where the second line uses the variance of the OLSE.

This, as shown in the answer that you link to (see also here), is independent of $$d:=(n-k)\frac{\hat{\sigma}^{2}}{\sigma^{2}}\sim\chi^{2}_{n-k},$$ where $$\hat{\sigma}^{2}=y'M_Xy/(n-k)$$ is the usual unbiased error variance estimate, with $$M_{X}=I-X(X'X)^{-1}X'$$ is the "residual maker matrix" from regressing on $$X$$.

So, as $$n'n$$ is a quadratic form in normals, $$\begin{eqnarray*} \frac{\overbrace{n^\prime n}^{\sim\chi^{2}_{q}}/q}{d/(n-k)}=\frac{(\hat{\beta}_{\text{ols}}-\beta)^\prime R\left\{R^\prime(X^\prime X)^{-1}R\right\}^{-1}R^\prime(\hat{\beta}_{\text{ols}}-\beta)/q}{\hat{\sigma}^{2}}\sim F_{q,n-k}. \end{eqnarray*}$$ In particular, under $$H_{0}:R^\prime\beta=r$$, this reduces to the statistic $$\begin{eqnarray} F=\frac{(R^\prime\hat{\beta}_{\text{ols}}-r)^\prime\left\{R^\prime(X^\prime X)^{-1}R\right\}^{-1}(R^\prime\hat{\beta}_{\text{ols}}-r)/q}{\hat{\sigma}^{2}}\sim F_{q,n-k}. \end{eqnarray}$$

For illustration, consider the special case $$R^\prime=I$$, $$r=0$$, $$q=2$$, $$\hat{\sigma}^{2}=1$$ and $$X^\prime X=I$$. Then, $$\begin{eqnarray} F=\hat{\beta}_{\text{ols}}^\prime\hat{\beta}_{\text{ols}}/2=\frac{\hat{\beta}_{\text{ols},1}^2+\hat{\beta}_{\text{ols},2}^2}{2}, \end{eqnarray}$$ the squared Euclidean distance of the OLS estimate from the origin standardized by the number of elements - highlighting that, since $$\hat{\beta}_{\text{ols},2}^2$$ are squared standard normals and hence $$\chi^2_1$$, the $$F$$ distribution may be seen as an "average $$\chi^2$$ distribution.

In case you prefer a little simulation (which is of course not a proof!), in which the null is tested that none of the $$k$$ regressors matter - which they indeed do not, so that we simulate the null distribution.

We see very good agreement between the theoretical density and the histogram of the Monte Carlo test statistics.

library(lmtest)
n <- 100
reps <- 20000
sloperegs <- 5 # number of slope regressors, q or k-1 (minus the
# constant) in the above notation
critical.value <- qf(p = .95, df1 = sloperegs, df2 = n-sloperegs-1)
# for the null that none of the slope regressors matter

Fstat <- rep(NA,reps)
for (i in 1:reps){
y <- rnorm(n)
X <- matrix(rnorm(n*sloperegs), ncol=sloperegs)
reg <- lm(y~X)
Fstat[i] <- waldtest(reg, test="F")F[2] } mean(Fstat>critical.value) # very close to 0.05 hist(Fstat, breaks = 60, col="lightblue", freq = F, xlim=c(0,4)) x <- seq(0,6,by=.1) lines(x, df(x, df1 = sloperegs, df2 = n-sloperegs-1), lwd=2, col="purple")  To see that the versions of the test statistics in the question and the answer are indeed equivalent, note that the null corresponds to the restrictions $$R'=[0\;\;I]$$ and $$r=0$$. Let $$X=[X_1\;\;X_2]$$ be partitioned according to which coefficients are restricted to be zero under the null (in your case, all but the constant, but the derivation to follow is general). Also, let $$\hat{\beta}_{\text{ols}}=(\hat{\beta}_{\text{ols},1}^\prime,\hat{\beta}_{\text{ols},2}')'$$ be the suitably partitioned OLS estimate. Then, $$R'\hat{\beta}_{\text{ols}}=\hat{\beta}_{\text{ols},2}$$ and $$R^\prime(X^\prime X)^{-1}R\equiv\tilde D,$$ the lower right block of \begin{align*} (X^TX)^{-1}&=\left( \begin{array} {c,c} X_1'X_1&X_1'X_2 \\ X_2'X_1&X_2'X_2\end{array} \right)^{-1}\\&\equiv\left( \begin{array} {c,c} \tilde A&\tilde B \\ \tilde C&\tilde D\end{array} \right) \end{align*} Now, use results for partitioned inverses to obtain $$\tilde D=(X_2'X_2-X_2'X_1(X_1'X_1)^{-1}X_1'X_2)^{-1}=(X_2'M_{X_1}X_2)^{-1}$$ where $$M_{X_1}=I-X_1(X_1'X_1)^{-1}X_1'$$. Thus, the numerator of the $$F$$ statistic becomes (without the division by $$q$$) $$F_{num}=\hat{\beta}_{\text{ols},2}'(X_2'M_{X_1}X_2)\hat{\beta}_{\text{ols},2}$$ Next, recall that by the Frisch-Waugh-Lovell theorem we may write $$\hat{\beta}_{\text{ols},2}=(X_2'M_{X_1}X_2)^{-1}X_2'M_{X_1}y$$ so that \begin{align*} F_{num}&=y'M_{X_1}X_2(X_2'M_{X_1}X_2)^{-1}(X_2'M_{X_1}X_2)(X_2'M_{X_1}X_2)^{-1}X_2'M_{X_1}y\\ &=y'M_{X_1}X_2(X_2'M_{X_1}X_2)^{-1}X_2'M_{X_1}y \end{align*} It remains to show that this numerator is identical to $$\text{RSSR}-\text{USSR}$$, the difference in restricted and unrestricted sum of squared residuals. Here, $$\text{RSSR}=y'M_{X_1}y$$ is the residual sum of squares from regressing $$y$$ on $$X_1$$, i.e., with $$H_0$$ imposed. In your special case, this is just $$TSS=\sum_i(y_i-\bar y)^2$$, the residuals of a regression on a constant. Again using FWL (which also shows that the residuals of the two approaches are identical), we can write $$\text{USSR}$$ (SSR in your notation) as the SSR of the regression $$M_{X_1}y\quad\text{on}\quad M_{X_1}X_2$$ That is, $$\begin{eqnarray*} \text{USSR}&=&y'M_{X_1}'M_{M_{X_1}X_2}M_{X_1}y\\ &=&y'M_{X_1}'(I-P_{M_{X_1}X_2})M_{X_1}y\\ &=&y'M_{X_1}y-y'M_{X_1}M_{X_1}X_2((M_{X_1}X_2)'M_{X_1}X_2)^{-1}(M_{X_1}X_2)'M_{X_1}y\\ &=&y'M_{X_1}y-y'M_{X_1}X_2(X_2'M_{X_1}X_2)^{-1}X_2'M_{X_1}y \end{eqnarray*}$$ Thus, $$\begin{eqnarray*} \text{RSSR}-\text{USSR}&=&y'M_{X_1}y-(y'M_{X_1}y-y'M_{X_1}X_2(X_2'M_{X_1}X_2)^{-1}X_2'M_{X_1}y)\\ &=&y'M_{X_1}X_2(X_2'M_{X_1}X_2)^{-1}X_2'M_{X_1}y \end{eqnarray*}$$ @ChristophHanck has provided a very comprehensive answer, here I will add a sketch of proof on the special case OP mentioned. Hopefully it's also easier to follow for beginners. A random variableY\sim F_{d_1,d_2}$if $$Y=\frac{X_1/d_1}{X_2/d_2},$$ where$X_1\sim\chi^2_{d_1}$and$X_2\sim\chi^2_{d_2}$are independent. Thus, to show that the$F$-statistic has$F$-distribution, we may as well show that$c\text{ESS}\sim\chi^2_{p-1}$and$c\text{RSS}\sim\chi^2_{n-p}$for some constant$c$, and that they are independent. In OLS model we write $$y=X\beta+\varepsilon,$$ where$X$is a$n\times p$matrix, and ideally$\varepsilon\sim N_n(\mathbf{0}, \sigma^2I)$. For convenience we introduce the hat matrix$H=X(X^TX)^{-1}X^{T}$(note$\hat{y}=Hy$), and the residual maker$M=I-H$. Important properties of$H$and$M$are that they are both symmetric and idempotent. In addition, we have$\operatorname{tr}(H)=p$and$HX=X$, these will come in handy later. Let us denote the matrix of all ones as$J$, the sum of squares can then be expressed with quadratic forms: $$\text{TSS}=y^T\left(I-\frac{1}{n}J\right)y,\quad\text{RSS}=y^TMy,\quad\text{ESS}=y^T\left(H-\frac{1}{n}J\right)y.$$ Note that$M+(H-J/n)+J/n=I$. One can verify that$J/n$is idempotent and$\operatorname{rank}(M)+\operatorname{rank}(H-J/n)+\operatorname{rank}(J/n)=n$. It follows from this then that$H-J/n$is also idempotent and$M(H-J/n)=0$. We can now set out to show that$F$-statistic has$F$-distribution (search Cochran's theorem for more). Here we need two facts: 1. Let$x\sim N_n(\mu,\Sigma)$. Suppose$A$is symmetric with rank$r$and$A\Sigma$is idempotent, then$x^TAx\sim\chi^2_r(\mu^TA\mu/2)$, i.e. non-central$\chi^2$with d.f.$r$and non-centrality$\mu^TA\mu/2$. This is a special case of Baldessari's result, a proof can also be found here. 2. Let$x\sim N_n(\mu,\Sigma)$. If$A\Sigma B=0$, then$x^TAx$and$x^TBx$are independent. This is known as Craig's theorem. Since$y\sim N_n(X\beta,\sigma^2I)$, we have $$\frac{\text{ESS}}{\sigma^2}=\left(\frac{y}{\sigma}\right)^T\left(H-\frac{1}{n}J\right)\frac{y}{\sigma}\sim\chi^2_{p-1}\left((X\beta)^T\left(H-\frac{J}{n}\right)X\beta\right).$$ However, under null hypothesis$\beta=\mathbf{0}$, so really$\text{ESS}/\sigma^2\sim\chi^2_{p-1}$. On the other hand, note that$y^TMy=\varepsilon^TM\varepsilon$since$HX=X$. Therefore$\text{RSS}/\sigma^2\sim\chi^2_{n-p}$. Since$M(H-J/n)=0$,$\text{ESS}/\sigma^2$and$\text{RSS}/\sigma^2$are also independent. It immediately follows then $$F = \frac{(\text{TSS}-\text{RSS})/(p-1)}{\text{RSS}/(n-p)}=\frac{\dfrac{\text{ESS}}{\sigma^2}/(p-1)}{\dfrac{\text{RSS}}{\sigma^2}/(n-p)}\sim F_{p-1,n-p}.$$ • shouldn't the null hypothesis be$\beta = [1 \space\space \mathbf{0} ]^T\$ ? Commented Sep 12, 2023 at 14:10
• @abhishek Yes, you are correct. Commented Mar 15 at 11:19

Christoph's answer derived the null distribution of $$F$$ by starting with analyzing the distribution of the OLS estimate $$\hat{\beta}_{\text{ols}}$$. Below I want to introduce an alternative approach which does not need to explicitly writing down the OLS estimate in $$F$$. To save some typing, without loss of generality, assume the variance of error $$\sigma^2 = 1$$ subsequently.

The key idea of this approach is by viewing the numerator (before scaled by $$p$$) $$TSS - RSS$$ of the $$F$$ statistic as $$RSS(M_{\text{reduced}}) - RSS(M_{\text{full}})$$, where $$M_{\text{full}}$$ is the "full model" (or the model without any constraints imposed on $$\beta$$) $$y = X\beta + \varepsilon$$, and $$M_{\text{reduced}}$$ is the "reduced model" as designated by the null hypothesis $$H_0: \beta_1 = \cdots = \beta_p = 0$$, that is, $$y = \beta_0e + \varepsilon$$ (where $$e$$ is an $$n$$-length vector of all ones). Also note that for any model $$M: \tilde{y} = \tilde{X}\beta + \varepsilon$$ such that $$\tilde{X} \in \mathbb{R}^{n \times q}$$ spans a subspace of the column space spanned by the full design matrix $$X$$ and $$\tilde{y}$$ is a translation of the original response vector $$y$$, its residual sum of squares $$RSS(M)$$ always equals to \begin{align*} RSS(M) = \tilde{y}^\top(I - H_{\tilde{X}})\tilde{y} = \varepsilon^\top(I - H_{\tilde{X}})\varepsilon, \tag{1}\label{1} \end{align*} where $$H_{\tilde{X}} = \tilde{X}(\tilde{X}^\top\tilde{X})^{-1}\tilde{X}^\top$$ is the well-known "hat-matrix" based on $$\tilde{X}$$. Note that the last equality in $$\eqref{1}$$ holds because $$\tilde{y} = \tilde{X}\beta + \varepsilon$$ and $$(I - H_{\tilde{X}})\tilde{X}\beta = \tilde{X}\beta - \tilde{X}\beta = 0$$.

By $$\eqref{1}$$, we have \begin{align*} & TSS - RSS = RSS(M_{\text{reduced}}) - RSS(M_{\text{full}}) \\ =& \varepsilon^\top(I - n^{-1}ee^\top)\varepsilon - \varepsilon^\top(I - X(X^\top X)^{-1}X^\top)\varepsilon \\ =:& \varepsilon^\top(H - n^{-1}ee^\top)\varepsilon \tag{2.1}\label{2.1}, \\ & RSS = RSS(M_{\text{full}}) = \varepsilon^\top(I - H)\varepsilon. \tag{2.2}\label{2.2} \end{align*} Next, verify that the matrix $$H - n^{-1}ee^\top$$ is symmetric and idempotent (by using the fact that $$He = e$$) with rank $$p - 1$$, and the matrix $$I - H$$ is symmetric and idempotent with rank $$n - p - 1$$. It follows by the theorem at the end of this answer that \begin{align*} & TSS - RSS = \varepsilon^\top(H - n^{-1}ee^\top)\varepsilon \sim \chi_{p - 1}^2, \\ & RSS = \varepsilon^\top(I - H)\varepsilon \sim \chi_{n - p - 1}^2. \tag{3}\label{3} \end{align*}

Finally, since $$(H - n^{-1}ee^\top)\varepsilon \sim N(0, H - n^{-1}ee^\top)$$ and $$(I - H)\varepsilon \sim N(0, I - H)$$ are independent (verify their correlation is $$0$$), we conclude that $$TSS - RSS$$ and $$RSS$$ are independent, hence $$\eqref{3}$$ and the definition of $$F$$-distribution immediately imply that \begin{align*} F = \frac{\frac{TSS - RSS}{p - 1}}{\frac{RSS}{n - p - 1}} \sim F_{p - 1, n - p - 1}. \end{align*}

This completes the proof.

The key idea employed by the above proof may be outlined as follows:

1. Simplify the full model $$M_{\text{full}}$$ to a reduced-form $$M_{\text{reduced}}$$ by "inserting" a solution $$\beta^*$$ to the linear system as specified by the null hypothesis $$H_0$$ into $$y = X\beta + \varepsilon$$.
2. Compute $$RSS(M_{\text{full}})$$ and $$RSS(M_{\text{reduced}})$$ respectively, which enables us to express the numerator and the denominator of an $$F$$-statistic in terms of proper quadratic forms. The degrees of freedom of this $$F$$-statistic are then easily determined by identifying ranks of these quadratic forms.

As a further illustration, let's find the null distribution of the $$F$$-statistic for testing the general linear hypothesis $$H_0: A\beta = \gamma$$ in model $$y = X\beta + \varepsilon$$. Here $$A \in \mathbb{R}^{k \times (p + 1)}$$ has rank $$r \leq \min(k, p + 1)$$, $$\gamma \in \mathbb{R}^k$$. By writing $$A\beta = \gamma$$, we tacitly assume that this linear system has at least one specific solution $$\beta^*$$ (otherwise $$H_0$$ would be trivially rejected).

To carry out step 1, note that by a structural theorem of the non-homogeneous linear system, a general solution to $$A\beta = \gamma$$ can be written as $$\beta = \beta^* + B\delta$$, where $$B \in \mathbb{R}^{(p + 1) \times (p + 1 - r)}$$ is of full column rank (specifically, $$B$$ is formed by $$p + 1 - r$$ basic solutions to the homogeneous linear system $$A\beta = 0$$), $$\delta \in \mathbb{R}^{p + 1 - r}$$ may be deemed as "new" free parameter of the reduced model. Inserting this $$\beta$$ into $$y = X\beta + \varepsilon$$, we have \begin{align*} y = X\beta^* + XB\delta + \varepsilon, \end{align*} or \begin{align*} \tilde{y} = \tilde{X}\delta + \varepsilon, \tag{4}\label{4} \end{align*} where $$\tilde{y} = y - X\beta^*, \tilde{X} = XB$$. Note that this reduced model $$\eqref{4}$$ itself can be deemed as a "full model" with a reduced free parameter $$\delta$$.

To carry out step 2, note that $$RSS(M_{\text{full}})$$ is identical to $$\eqref{2.2}$$ (which is independent of the form of $$H_0$$). In addition, by $$\eqref{1}$$, we have \begin{align*} RSS(M_{\text{reduced}}) = \varepsilon^\top(I - H_{\tilde{X}})\varepsilon, \end{align*} where $$H_{\tilde{X}} = \tilde{X}(\tilde{X}^\top\tilde{X})^{-1}\tilde{X}^\top = XB(B^\top X^\top XB)^{-1}B^\top X^\top$$. It then follows that \begin{align*} RSS(M_{\text{reduced}}) - RSS(M_{\text{full}}) = \varepsilon^\top(X(X^\top X)^{-1}X^\top - XB(B^\top X^\top XB)^{-1}B^\top X^\top)\varepsilon =: \varepsilon^\top C \varepsilon. \end{align*} It is easy to verify that $$C$$ is symmetric and idempotent, whence \begin{align*} & \operatorname{rank}(C) = \operatorname{tr}(C) \\ =& \operatorname{tr}(X(X^\top X)^{-1}X^\top) - \operatorname{tr}(XB(B^\top X^\top XB)^{-1}B^\top X^\top) \\ =& \operatorname{tr}(X^\top X(X^\top X)^{-1}) - \operatorname{tr}(B^\top X^\top XB(B^\top X^\top XB)^{-1}) \\ =& \operatorname{tr}(I_{(p + 1)}) - \operatorname{tr}(I_{(p + 1 - r)}) \\ =& p + 1 - (p + 1 - r) = r. \end{align*} Therefore, under $$H_0$$, $$RSS(M_{\text{reduced}}) - RSS(M_{\text{full}}) = \varepsilon^\top C \varepsilon \sim \chi_r^2$$.

Furthermore, in view of \begin{align*} & \operatorname{Cov}(C\varepsilon, (I - H)\varepsilon) = C(I - H) \\ =& (H - XB(B^\top X^\top XB)^{-1}B^\top X^\top)(I - H) \\ =& H - XB(B^\top X^\top XB)^{-1}B^\top X^\top - H + XB(B^\top X^\top XB)^{-1}B^\top X^\top X(X^\top X)^{-1}X^\top \\ =& H - XB(B^\top X^\top XB)^{-1}B^\top X^\top - H + XB(B^\top X^\top XB)^{-1}B^\top X^\top = 0, \end{align*} $$RSS(M_{\text{reduced}}) - RSS(M_{\text{full}})$$ and $$RSS(M_{\text{full}})$$ are independent.

In conclusion, the null distribution of $$F$$ in this case is $$F_{r, n - p - 1}$$.

Reference Notes

I learned this approach from a Chinese linear model monograph 近代回归分析 ("Modern Regression Analysis") by Xiru Chen, Songgui Wang. The first author Chen is well recognized as one of the greatest Chinese mathematical statisticians.

I believe the treatment presented in Michael J. Wichura's text The Coordinate-Free Approach to Linear Models (check Theorem 3.8.2 and Chapter 8) is of the same spirit, but is of more geometric flavor (for example, the quadratic form $$\varepsilon^\top (I - H)\varepsilon$$ is always presented as $$\|P_{\operatorname{col}(X)^\perp}\varepsilon\|^2$$, where "$$P$$" stands for the orthogonal projection operator). Although his unique writing style looks intimidating at the first reading, the whole book is full of insights and eye-opening (e.g., many standard lengthy algebraic proofs are shortened significantly by his geometric argument), and it will surely bring your understanding of linear models to the next level.

• Nice and comprehensive (as usual). If I am not wrong, the same approach has been used by Searle in his text (and almost similar path but bit handwavy in Scheffé) and in more modern text by Christensen. Also, on a side note, could you be able to provide a link to the monograph by Xiru Chen? Embarrassingly I was not aware about any such work. Commented Mar 16 at 3:35
• @User1865345 Thanks for providing more references on this important topic. For Chen's book's link, I couldn't find (an English) one. This google book's link seems not providing much useful information. This book probably doesn't have an English translation. Commented Mar 16 at 3:43
• @User1865345 Chen wrote a sequence of textbooks/monographs (notably about 10 books that I am aware of), covering most classical areas in theoretical statistics. One that is very unique IMO is called "Advanced Mathematical Statistics" 高等数理统计学, in which more than half of the volume are devoted to (very challenging) problems (I have never seen a book has such breadth and depth of exercises). If you are interested, I could send pdfs to you by email (then maybe with the help of Chatgpt, you can dabble in reading it :)). Commented Mar 16 at 3:58
• I will not continue this discussion here, as it is getting off topic germane to the post. If you have time, you can visit Ten Fold and I more or less drop by regularly. Commented Mar 16 at 4:21
• @User1865345 Surely I will. Commented Mar 16 at 4:31