多元回归研究的是因变量 Y Y Y与 m m m个自变量 x 1 , ⋯   , x m x_1,\cdots,x_m x1​,⋯,xm​之间的相互依赖关系的，这里假定 Y Y Y是随机变量，而 x 1 , ⋯   , x m x_1,\cdots,x_m x1​,⋯,xm​是一般变量（视为非随机的）， Y Y Y不能由 x 1 , ⋯   , x m x_1,\cdots,x_m x1​,⋯,xm​完全决定，但是可以用 x 1 , ⋯   , x m x_1,\cdots,x_m x1​,⋯,xm​的函数进行估计。而多元线性回归中的，指的是关于“参数”的线性，我们一般可以把数据写成这样的线性模型： { y 1 = β 0 + β 1 x 11 + ⋯ + β m x 1 m + ε 1 , ⋯ ⋯ ⋯ ⋯ ⋯ ⋯ y n = β 0 + β 1 x n 1 + ⋯ + β m x n m + ε n . \left\{ \begin{array}l y_1=\beta_0+\beta_1x_{11}+\cdots+\beta_mx_{1m}+\varepsilon_1, \\ \cdots \cdots \cdots \cdots \cdots \cdots \\ y_n =\beta_0+\beta_1x_{n1}+\cdots+\beta_mx_{nm}+\varepsilon_n. \end{array} \right. ⎩⎨⎧​y1​=β0​+β1​x11​+⋯+βm​x1m​+ε1​,⋯⋯⋯⋯⋯⋯yn​=β0​+β1​xn1​+⋯+βm​xnm​+εn​.​ 这就是经典多元线性回归模型，这里 ε t \varepsilon_t εt​是偏差，可以将其视为一个与变量无关的白噪声，满足 E ( ε t ) = 0 , D ( ε t ) = σ 2 , C o v ( ε i , ε j ) = 0 ( i ≠ j ) {\rm E}(\varepsilon_t)=0,{\rm D}(\varepsilon_t)=\sigma^2,{\rm Cov}(\varepsilon_i,\varepsilon_j)=0(i\ne j) E(εt​)=0,D(εt​)=σ2,Cov(εi​,εj​)=0(i​=j)，还可以进一步假定其为正态白噪声，即 ε t ∼ N ( 0 , σ 2 ) \varepsilon_t\sim N(0,\sigma^2) εt​∼N(0,σ2)相互独立。

大型的式子都可以用向量、矩阵形式表达，多元线性回归模型也不例外，这里进行以下定义： Y = [ y 1 y 2 ⋮ y n ] n , β = [ β 0 β 1 ⋮ β m ] m + 1 , ε = [ ε 1 ε 2 ⋮ ε n ] n , C = ( 1 m ∣ X ) = [ 1 x 11 x 12 ⋯ x 1 m 1 x 21 x 22 ⋯ x 2 m ⋮ ⋮ ⋮ ⋮ 1 x n 1 x n 2 ⋯ x n m ] n × m , Y=\begin{bmatrix} y_1 \\ y_2 \\ \vdots \\ y_n \end{bmatrix}_n,\beta=\begin{bmatrix} \beta_0 \\ \beta_1 \\ \vdots \\ \beta_m \end{bmatrix}_{m+1},\varepsilon=\begin{bmatrix} \varepsilon_1 \\ \varepsilon_2 \\ \vdots \\ \varepsilon_n \end{bmatrix}_n,C=(\boldsymbol 1_m|X)=\begin{bmatrix} 1 & x_{11} & x_{12} & \cdots & x_{1m} \\ 1 & x_{21} & x_{22} & \cdots & x_{2m} \\ \vdots & \vdots & \vdots & & \vdots \\ 1 & x_{n1} & x_{n2} & \cdots & x_{nm} \end{bmatrix}_{n\times m}, Y=⎣⎢⎢⎢⎡​y1​y2​⋮yn​​⎦⎥⎥⎥⎤​n​,β=⎣⎢⎢⎢⎡​β0​β1​⋮βm​​⎦⎥⎥⎥⎤​m+1​,ε=⎣⎢⎢⎢⎡​ε1​ε2​⋮εn​​⎦⎥⎥⎥⎤​n​,C=(1m​∣X)=⎣⎢⎢⎢⎡​11⋮1​x11​x21​⋮xn1​​x12​x22​⋮xn2​​⋯⋯⋯​x1m​x2m​⋮xnm​​⎦⎥⎥⎥⎤​n×m​, 这样就有 Y = C β + ε Y=C\beta +\varepsilon Y=Cβ+ε，其中 ε ∼ N n ( 0 , σ 2 I n ) \varepsilon\sim N_n(0,\sigma^2I_n) ε∼Nn​(0,σ2In​)。要能够估计出 m + 1 m+1 m+1个参数 β 0 , ⋯   , β m \beta_0,\cdots,\beta_m β0​,⋯,βm​，就至少要有 m + 1 m+1 m+1个样本，也就是 n > m n>m n>m；同时我们要求各个自变量之间是不相关的，因为如果存在相关性，则自变量可以由其他自变量表出，所需要的回归系数就可以减少，所以我们的要求是 r a n k ( C ) = m + 1 {\rm rank} (C)=m+1 rank(C)=m+1。

在给定模型后，我们的目的是估计回归系数 β \beta β，在多元线性回归中常用的方法是最小二乘法。最小二乘法的思想就是让偏差平方和达到最小，偏差平方和即 Q ( β ) = ∑ t = 1 n e t 2 = ∑ t = 1 n [ y t − ( β 0 + β 1 x t 1 + ⋯ + β m x m t ) ] 2 = ( Y − C β ) ′ ( Y − C β ) . Q(\beta)=\sum_{t=1}^n e_t^2=\sum_{t=1}^n [y_t-(\beta_0+\beta_1x_{t1}+\cdots+\beta_mx_{mt})]^2=(Y-C\beta)'(Y-C\beta). Q(β)=t=1∑n​et2​=t=1∑n​[yt​−(β0​+β1​xt1​+⋯+βm​xmt​)]2=(Y−Cβ)′(Y−Cβ).

如果令 Y t Y_t Yt​为 Y Y Y的第 t t t行， C t C_t Ct​为 C C C的第 t t t行，则 y t − ( β 0 + β 1 x t 1 + ⋯ + β m x m t ) = Y t − C t β , [ y t − ( β 0 + β 1 x t 1 + ⋯ + β m x m t ) ] 2 = ( Y t − C t β ) ′ ( Y t − C t β ) . y_t-(\beta_0+\beta_1x_{t1}+\cdots+\beta_mx_{mt})=Y_t-C_t\beta,\\ [y_t-(\beta_0+\beta_1x_{t1}+\cdots+\beta_mx_{mt})]^2=(Y_t-C_t\beta)'(Y_t-C_t\beta). yt​−(β0​+β1​xt1​+⋯+βm​xmt​)=Yt​−Ct​β,[yt​−(β0​+β1​xt1​+⋯+βm​xmt​)]2=(Yt​−Ct​β)′(Yt​−Ct​β). 所以 ( Y − C β ) ′ ( Y − C β ) = [ ( Y 1 − C 1 β ) ′ ( Y 2 − C 2 β ) ′ ⋮ ( Y n − C n β ) ′ ] [ ( Y 1 − C 1 β ) ( Y 2 − C 2 β ) ⋯ ( Y n − C n β ) ] = ∑ t = 1 n ( Y t − C t β ) ′ ( Y t − C t β ) = ∑ t = 1 n e t 2 . \begin{aligned} & (Y-C\beta)'(Y-C\beta) \\ =& \begin{bmatrix} (Y_1-C_1\beta)' \\ (Y_2-C_2\beta)' \\ \vdots \\ (Y_n-C_n\beta)' \end{bmatrix}\begin{bmatrix} (Y_1-C_1\beta) & (Y_2-C_2\beta) & \cdots & (Y_n-C_n\beta) \end{bmatrix}\\ =&\sum_{t=1}^n (Y_t-C_t\beta)'(Y_t-C_t\beta) \\ =&\sum_{t=1}^n e_t^2. \end{aligned} ===​(Y−Cβ)′(Y−Cβ)⎣⎢⎢⎢⎡​(Y1​−C1​β)′(Y2​−C2​β)′⋮(Yn​−Cn​β)′​⎦⎥⎥⎥⎤​[(Y1​−C1​β)​(Y2​−C2​β)​⋯​(Yn​−Cn​β)​]t=1∑n​(Yt​−Ct​β)′(Yt​−Ct​β)t=1∑n​et2​.​

要使得 Q ( β ) Q(\beta) Q(β)最小，就有 ∂ Q ( β ) ∂ β = ∂ ( Y − C β ) ′ ( Y − C β ) ∂ β = ∂ ( Y ′ Y ) ∂ β − ∂ ( β ′ C ′ Y ) ∂ β − ∂ ( Y ′ C β ) ∂ β + ∂ ( β ′ C ′ C β ) ∂ β = 0 m + 1 − C ′ Y − C ′ Y + 2 C ′ C β = 2 ( C ′ C β − C ′ Y ) \begin{aligned} \frac{\partial Q(\beta)}{\partial \beta}=& \frac{\partial (Y-C\beta)'(Y-C\beta)}{\partial \beta} \\ =&\frac{\partial(Y'Y)}{\partial \beta}-\frac{\partial(\beta'C'Y)}{\partial \beta}-\frac{\partial(Y'C\beta)}{\partial \beta}+\frac{\partial(\beta'C'C\beta)}{\partial \beta}\\ =&\boldsymbol 0_{m+1}-C'Y-C'Y+2C'C\beta\\ =&2(C'C\beta-C'Y) \end{aligned} ∂β∂Q(β)​====​∂β∂(Y−Cβ)′(Y−Cβ)​∂β∂(Y′Y)​−∂β∂(β′C′Y)​−∂β∂(Y′Cβ)​+∂β∂(β′C′Cβ)​0m+1​−C′Y−C′Y+2C′Cβ2(C′Cβ−C′Y)​

这里用到的矩阵微商公式有： ∂ A ∂ β = 0 , ∂ β ′ y ∂ β = y , ∂ y ′ β ∂ β = y 或 ∂ y ′ A β ∂ β = A ′ y , ∂ β ′ A β ∂ β = ( A + A ′ ) β . \frac{\partial A}{\partial \beta}=0, \\ \frac{\partial \beta'y}{\partial \beta}=y, \\ \frac{\partial y'\beta}{\partial\beta}=y或\frac{\partial y'A\beta}{\partial \beta}=A'y,\\ \frac{\partial \beta'A\beta}{\partial \beta}=(A+A')\beta. ∂β∂A​=0,∂β∂β′y​=y,∂β∂y′β​=y或∂β∂y′Aβ​=A′y,∂β∂β′Aβ​=(A+A′)β.

所以让导数为0，就得到 β ^ = ( C ′ C ) − 1 C ′ Y . \hat \beta=(C'C)^{-1}C'Y. β^​=(C′C)−1C′Y. 在刚才的推导中，我们得到了最小二乘法估计的参数值 β ^ = ( C ′ C ) − 1 C ′ Y \hat\beta=(C'C)^{-1}C'Y β^​=(C′C)−1C′Y，结合我们的多元线性回归模型 Y = C β + ε Y=C\beta+\varepsilon Y=Cβ+ε，在已知 X X X时 Y Y Y的预测值就应该是 Y ^ = C β ^ = C ( C ′ C ) − 1 C ′ Y = d H Y . \hat Y=C\hat \beta=C(C'C)^{-1}C' Y\stackrel {\rm d}=HY. Y^=Cβ^​=C(C′C)−1C′Y=dHY. 这里的 H = C ( C ′ C ) − 1 C ′ H=C(C'C)^{-1}C' H=C(C′C)−1C′被称为帽子矩阵，可以看到 H ′ H = C ( C ′ C ) − 1 C ′ ⋅ C ( C ′ C ) − 1 C ′ = H 2 = C ( C ′ C ) − 1 C ′ = H . H'H=C(C'C)^{-1}C'\cdot C(C'C)^{-1}C'=H^2=C(C'C)^{-1}C'=H. H′H=C(C′C)−1C′⋅C(C′C)−1C′=H2=C(C′C)−1C′=H. 并且，此时的偏差为 e = ε ^ = Y − Y ^ = ( I − H ) Y e=\hat\varepsilon=Y-\hat Y=(I-H)Y e=ε^=Y−Y^=(I−H)Y，残差平方和就是 Q ( β ^ ) = ε ^ ′ ε ^ = Y ′ ( I − H ) ′ ( I − H ) Y = Y ′ ( I − H ) Y = Y ′ Y − Y ′ H Y = Y ′ Y − Y ′ C β ^ . \begin{aligned} Q(\hat\beta)=&\hat\varepsilon'\hat\varepsilon \\ =&Y'(I-H)'(I-H)Y \\ =&Y'(I-H)Y\\ =&Y'Y-Y'HY\\ =&Y'Y-Y'C\hat \beta. \end{aligned} Q(β^​)=====​ε^′ε^Y′(I−H)′(I−H)YY′(I−H)YY′Y−Y′HYY′Y−Y′Cβ^​.​

最小二乘估计量 β ^ = ( C ′ C ) − 1 C ′ Y \hat \beta=(C'C)^{-1}C'Y β^​=(C′C)−1C′Y是 β \beta β的最小方差线性无偏估计(BLUE)，对这个词做一些说明。

首先，很明显 β ^ \hat \beta β^​是一个随机变量（因为 Y Y Y被假定为随机向量， X X X被视为普通向量，所以 C C C应该被视为常数阵， Y Y Y视为随机向量），而线性性，指的就是 β ^ \hat \beta β^​可以被 Y 1 , ⋯   , Y n Y_1,\cdots,Y_n Y1​,⋯,Yn​所线性表示，这里的 ( C ′ C ) − 1 C ′ (C'C)^{-1}C' (C′C)−1C′就是其线性表示系数。线性性的好处是，一些对于随机变量数字特征的对于独立随机变量的线性情形是有计算公式的，这样就可以很好地写出 β ^ \hat \beta β^​的表现形式及其数字特征。

其次，无偏性指的是 E β ^ = β {\rm E}\hat \beta=\beta Eβ^​=β，注意这里的 β \beta β是客观存在的，是不随抽样改变的常数，而 β ^ \hat \beta β^​是随着我们抽取样本进行观测，根据样本的观测值算出的统计量，它是一个具有两重性的统计量。而 Y = C β + ε , E ( Y ) = C β , E β ^ = ( C ′ C ) − 1 C ′ E ( Y ) = ( C ′ C ) − 1 ( C ′ C ) β = β . D β ^ = E [ ( β ^ − β ) ( β ^ − β ) ′ ] = E [ ( C ′ C ) − 1 C ′ μ μ ′ C ( C ′ C ) − 1 ] = ( C ′ C ) − 1 C ′ E μ μ ′ C ( C ′ C ) − 1 = σ 2 ( C ′ C ) − 1 . Y=C\beta+\varepsilon,\quad {\rm E}(Y)=C\beta, \\ {\rm E}\hat \beta=(C'C)^{-1}C'{\rm E}(Y)=(C'C)^{-1}(C'C)\beta=\beta. \\ \begin{aligned} {\rm D}\hat \beta=&{\rm E}[(\hat\beta-\beta)(\hat\beta-\beta)'] \\ =&{\rm E}[(C'C)^{-1}C'\mu\mu'C(C'C)^{-1}]\\ =&(C'C)^{-1}C'{\rm E}\mu\mu' C(C'C)^{-1}\\ =&\sigma^2(C'C)^{-1}. \end{aligned} Y=Cβ+ε,E(Y)=Cβ,Eβ^​=(C′C)−1C′E(Y)=(C′C)−1(C′C)β=β.Dβ^​====​E[(β^​−β)(β^​−β)′]E[(C′C)−1C′μμ′C(C′C)−1](C′C)−1C′Eμμ′C(C′C)−1σ2(C′C)−1.​ 最后，最小方差性，指的是在所有的线性无偏估计中， β ^ \hat \beta β^​是方差最小的一个。假设还有其他的线性无偏估计量 β ∗ \beta^* β∗，记作 [ ( C ′ C ) − 1 C ′ + D ] Y = C ∗ Y = β ^ + D Y [(C'C)^{-1}C'+D]Y=C^*Y=\hat \beta+DY [(C′C)−1C′+D]Y=C∗Y=β^​+DY，为了满足无偏性，有 E β ∗ = E [ C ∗ ( X β + μ ) ] = E ( C ∗ X ) β = β , C ∗ X = I . {\rm E}\beta^*={\rm E}[C^*(X\beta+\mu)]={\rm E}(C^*X)\beta=\beta,\quad C^*X=I. Eβ∗=E[C∗(Xβ+μ)]=E(C∗X)β=β,C∗X=I. 而 C ∗ X = ( X ′ X ) − 1 X ′ X + D X = I + D X = I C^*X=(X'X)^{-1}X'X+DX=I+DX=I C∗X=(X′X)−1X′X+DX=I+DX=I，所以 D X = 0 DX=0 DX=0。于是 C o v ( β ∗ ∣ X ) = E [ ( C ∗ Y − β ) ( C ∗ Y − β ) ′ ∣ X ] = [ E ( C ∗ μ ) ( C ∗ μ ) ′ ∣ X ] = E [ ( C ∗ μ μ ′ C ∗ ′ ) ∣ X ] = [ ( X ′ X ) − 1 X ′ + D ] E ( μ μ ′ ∣ X ) [ X ( X ′ X ) − 1 + D ′ ] = σ 2 [ ( X ′ X ) − 1 X ′ X ( X ′ X ) − 1 + ( X ′ X ) − 1 ( X D ) ′ + D X ( X ′ X ) − 1 + D D ′ ] = σ 2 ( X ′ X ) − 1 + σ 2 D D ′ . \begin{aligned} {\rm Cov}(\beta^*|X)=&{\rm E}[(C^*Y-\beta)(C^*Y-\beta)'|X]\\ =&[{\rm E}(C^*\mu)(C^*\mu)'|X]\\ =&{\rm E}[(C^*\mu\mu'{C^*}')|X]\\ =&[(X'X)^{-1}X'+D]{\rm E}(\mu\mu'|X)[X(X'X)^{-1}+D']\\ =&\sigma^2[(X'X)^{-1}X'X(X'X)^{-1}+(X'X)^{-1}(XD)'+DX(X'X)^{-1}+DD']\\ =&\sigma^2(X'X)^{-1}+\sigma^2DD'. \end{aligned} Cov(β∗∣X)======​E[(C∗Y−β)(C∗Y−β)′∣X][E(C∗μ)(C∗μ)′∣X]E[(C∗μμ′C∗′)∣X][(X′X)−1X′+D]E(μμ′∣X)[X(X′X)−1+D′]σ2[(X′X)−1X′X(X′X)−1+(X′X)−1(XD)′+DX(X′X)−1+DD′]σ2(X′X)−1+σ2DD′.​ 由于 D D D是正定的，所以 σ 2 D D ′ > 0 \sigma^2DD'>0 σ2DD′>0，所以 C o v ( β ∗ ∣ X ) > σ 2 ( X ′ X ) − 1 = C o v ( β ^ ∣ X ) {\rm Cov}(\beta^*|X)>\sigma^2(X'X)^{-1}={\rm Cov}(\hat \beta|X) Cov(β∗∣X)>σ2(X′X)−1=Cov(β^​∣X)，这就证明了 β ^ \hat \beta β^​的最小方差性。 这就说明 β ^ \hat \beta β^​是最小方差的。

因为 β ^ \hat \beta β^​是 β \beta β的最小方差线性无偏估计，所以对于 β \beta β的线性组合 α ′ β \alpha'\beta α′β，其最小方差线性无偏估计依然是 a ′ β ^ a'\hat\beta a′β^​。实际应用时，我们会得到一组数 x t 1 , ⋯   , x t m x_{t1},\cdots,x_{tm} xt1​,⋯,xtm​，要用它预报 y t y_t yt​，实际上就是对 β \beta β分量的线性组合，这也是我们用 C β ^ C\hat\beta Cβ^​作为 Y Y Y的预测值的理由。

接下来探讨 β ^ \hat \beta β^​作为一个随机变量时的分布。首先由线性性， β ^ \hat\beta β^​是 Y Y Y的分量的线性组合，而 Y Y Y的每一个分量中随机部分都是 ε \varepsilon ε分量，在我们假定 ε \varepsilon ε是正态白噪声的情况下， Y Y Y的各个分量应该是相互独立且服从正态分布的，即 Y ∼ N n ( C β , σ 2 I n ) Y\sim N_n(C\beta,\sigma^2I_n) Y∼Nn​(Cβ,σ2In​)，所以 E ( β ^ ) = ( C ′ C ) − 1 C ′ E ( Y ) = β , D ( β ^ ) = ( C ′ C ) − 1 C ′ D ( Y ) C ( C ′ C ) − 1 = σ 2 ( C ′ C ) − 1 . {\rm E}(\hat\beta)=(C'C)^{-1}C'{\rm E}(Y)=\beta, \\ {\rm D}(\hat\beta)=(C'C)^{-1}C'{\rm D}(Y)C(C'C)^{-1}=\sigma^2(C'C)^{-1}. E(β^​)=(C′C)−1C′E(Y)=β,D(β^​)=(C′C)−1C′D(Y)C(C′C)−1=σ2(C′C)−1. 即 β ^ ∼ N m + 1 ( β , σ 2 ( C ′ C ) − 1 ) \hat \beta\sim N_{m+1}(\beta,\sigma^2(C'C)^{-1}) β^​∼Nm+1​(β,σ2(C′C)−1)。

最后简要提一下白噪声方差 σ 2 \sigma^2 σ2的估计，由最大似然估计法， β \beta β的极大似然估计依然是 β ^ \hat \beta β^​，还可以得到 σ 2 \sigma^2 σ2的最大似然估计是 σ ^ 2 = ( Y − C β ^ ) ′ ( Y − C β ^ ) n = Q ( β ^ ) n . \hat \sigma^2=\frac{(Y-C\hat\beta)'(Y-C\hat\beta)}{n}=\frac{Q(\hat\beta)}{n}. σ^2=n(Y−Cβ^​)′(Y−Cβ^​)​=nQ(β^​)​. 但这个估计量不是无偏的，所以我们对其作修正，就得到了最小二乘估计 s 2 = ( Y − C β ^ ) ′ ( Y − C β ^ ) n − m − 1 , E ( s 2 ) = σ 2 . s^2=\frac{(Y-C\hat\beta)'(Y-C\hat\beta)}{n-m-1},\quad {\rm E}(s^2)=\sigma^2. s2=n−m−1(Y−Cβ^​)′(Y−Cβ^​)​,E(s2)=σ2.

多元线性回归模型指的是这样的模型： Y = C β + ε Y=C\beta+\varepsilon Y=Cβ+ε，其中 C = ( 1 n ∣ X ) C=(\boldsymbol 1_n|X) C=(1n​∣X)， ε t ∼ W N ( 0 , σ 2 ) \varepsilon_t\sim {\rm WN}(0,\sigma^2) εt​∼WN(0,σ2)。即使 ε t \varepsilon_t εt​不是正态白噪声，也至少应当满足零均值、同方差、序列不相关性。

使用最小二乘法求解多元线性回归模型，得到 β \beta β的最小二乘估计量是 β ^ = ( C ′ C ) − 1 C ′ Y \hat\beta=(C'C)^{-1}C'Y β^​=(C′C)−1C′Y，基于此估计量和数据 C C C估计 Y ^ = C β ^ = C ( C ′ C ) − 1 C ′ Y = H Y \hat Y=C\hat\beta=C(C'C)^{-1}C'Y=HY Y^=Cβ^​=C(C′C)−1C′Y=HY，将 H H H称为帽子矩阵。

β ^ \hat\beta β^​是 β \beta β的最佳线性无偏估计(BLUE)，即满足线性性、无偏性、最小方差性。且其分布是 β ^ ∼ N m + 1 ( β , σ 2 ( C ′ C ) − 1 ) . \hat\beta\sim N_{m+1}(\beta,\sigma^2(C'C)^{-1}). β^​∼Nm+1​(β,σ2(C′C)−1).

在获得 β ^ \hat\beta β^​后，残差平方和为 Q ( β ^ ) = ( Y − C β ^ ) ′ ( Y − C β ^ ) = Y ′ Y − Y ′ C β ^ Q(\hat\beta)=(Y-C\hat\beta)'(Y-C\hat\beta)=Y'Y-Y'C\hat\beta Q(β^​)=(Y−Cβ^​)′(Y−Cβ^​)=Y′Y−Y′Cβ^​，可以用残差平方和估计 σ 2 \sigma^2 σ2，其无偏估计为 s 2 = Q ( β ^ ) / ( n − m − 1 ) s^2=Q(\hat\beta)/(n-m-1) s2=Q(β^​)/(n−m−1)。