Let lowercase bold letters denote vectors. The linear model is:

$$ y_i = \mathbf{x}_i' \mathbf{b} + \epsilon_i $$

In matrix notation for all observations $i=1, \ldots, n$:

$$\mathbf{y} = X \mathbf{b} + \boldsymbol{\epsilon} $$

The OLS estimator for $\mathbf{b}$ is:

$$\hat{\mathbf{b}} = (X'X)^{-1}X'\mathbf{y} $$

Substituting:

$$\hat{\mathbf{b}} = (X'X)^{-1}X' \left( X \mathbf{b} + \boldsymbol{\epsilon}\right) $$

What's the variance of our estimator $\hat{\mathbf{b}}$? Take the variance (conditional on X)$:

\begin{align*} \operatorname{Var}\left( \hat{\mathbf{b}} \mid X\right) &= \operatorname{Var}\left( (X'X)^{-1}X' \left( X \mathbf{b} + \boldsymbol{\epsilon}\right) \mid X\right) \\ &= \operatorname{Var}\left( (X'X)^{-1}X' \boldsymbol{\epsilon} \mid X\right) \\ &= (X'X)^{-1}X' \operatorname{Var}\left( \boldsymbol{\epsilon} \mid X\right) X(X'X)^{-1} \end{align*}

We assumed that that $\epsilon_i$ were IID. Hence for all $i$ we have $\operatorname{Var}(\epsilon_i) = \sigma^2$, and for the whole vector $\boldsymbol{\epsilon})$ we have $\operatorname{Var}(\boldsymbol{\epsilon}) = \sigma^2 I$ where $I$ is the identity matrix.

Continuing: \begin{align*} \operatorname{Var}\left( \hat{\mathbf{b}} \mid X\right) &= (X'X)^{-1}X' \sigma^2 IX(X'X)^{-1} \\ &= \sigma^2 (X'X)^{-1}X'X(X'X)^{-1} \\ &= \sigma^2 (X'X)^{-1} \end{align*}

$\sigma^2$ is a scalar and can be moved wherever by the commutative property of multiplication. $(X'X)^{-1}X'X = I$ by definition of an inverse matrix.

You can estimate $\sigma^2$ using the residuals from your OLS regression. Let $e_i$ be the residual (as opposed to error-term) for observation $i$. That is:

$$ e_i = y_i - \mathbf{x}_i' \hat{\mathbf{b}} $$

Then the usual estimator for $\sigma^2$ is: $$\hat{\sigma}^2 = \frac{1}{n-k} \sum_i e_i^2$$ where $k$ is your number of regressors (including a constant).

Let lowercase bold letters denote vectors. The linear model is:

$$ y_i = \mathbf{x}_i' \mathbf{b} + \epsilon_i $$

In matrix notation for all observations $i=1, \ldots, n$:

$$\mathbf{y} = X \mathbf{b} + \boldsymbol{\epsilon} $$

The OLS estimator for $\mathbf{b}$ is:

$$\hat{\mathbf{b}} = (X'X)^{-1}X'\mathbf{y} $$

Substituting:

$$\hat{\mathbf{b}} = (X'X)^{-1}X' \left( X \mathbf{b} + \boldsymbol{\epsilon}\right) $$

What's the variance of our estimator $\hat{\mathbf{b}}$? Take the variance (conditional on X)$:

\begin{align*} \operatorname{Var}\left( \hat{\mathbf{b}} \mid X\right) &= \operatorname{Var}\left( (X'X)^{-1}X' \left( X \mathbf{b} + \boldsymbol{\epsilon}\right) \mid X\right) \\ &= \operatorname{Var}\left( (X'X)^{-1}X' \boldsymbol{\epsilon} \mid X\right) \\ &= (X'X)^{-1}X' \operatorname{Var}\left( \boldsymbol{\epsilon} \mid X\right) X(X'X)^{-1} \end{align*}

We assumed that that $\epsilon_i$ were IID. Hence for all $i$ we have $\operatorname{Var}(\epsilon_i) = \sigma^2$, and for the whole vector $\boldsymbol{\epsilon}$ we have $\operatorname{Var}(\boldsymbol{\epsilon}) = \sigma^2 I$ where $I$ is the identity matrix.

Continuing: \begin{align*} \operatorname{Var}\left( \hat{\mathbf{b}} \mid X\right) &= (X'X)^{-1}X' \sigma^2 IX(X'X)^{-1} \\ &= \sigma^2 (X'X)^{-1}X'X(X'X)^{-1} \\ &= \sigma^2 (X'X)^{-1} \end{align*}

$\sigma^2$ is a scalar and can be moved wherever by the commutative property of multiplication. $(X'X)^{-1}X'X = I$ by definition of an inverse matrix.

You can estimate $\sigma^2$ using the residuals from your OLS regression. Let $\hat{e}_i$ be the residual (as opposed to error-term $\epsilon_i$) for observation $i$. That is:

$$ \hat{e}_i = y_i - \mathbf{x}_i' \hat{\mathbf{b}} $$ Note that this is based upon estimate $\hat{\mathbf{b}}$ instead of true values $\mathbf{b}$. Then the usual estimator for $\sigma^2$ is: $$\hat{\sigma}^2 = \frac{1}{n-k} \sum_i \hat{e}_i^2$$ where $k$ is your number of regressors (including a constant).

Let lowercase bold letters denote vectors. The linear model is:

$$ y_i = \mathbf{x}_i' \mathbf{b} + \epsilon_i $$

In matrix notation for all observations $i=1, \ldots, n$:

$$\mathbf{y} = X \mathbf{b} + \boldsymbol{\epsilon} $$

The OLS estimator for $\mathbf{b}$ is:

$$\hat{\mathbf{b}} = (X'X)^{-1}X'\mathbf{y} $$

Substituting:

$$\hat{\mathbf{b}} = (X'X)^{-1}X' \left( X \mathbf{b} + \boldsymbol{\epsilon}\right) $$

Take the variance (conditional on X)$:

\begin{align*} \operatorname{Var}\left( \hat{\mathbf{b}} \mid X\right) &= \operatorname{Var}\left( (X'X)^{-1}X' \left( X \mathbf{b} + \boldsymbol{\epsilon}\right) \mid X\right) \\ &= \operatorname{Var}\left( (X'X)^{-1}X' \boldsymbol{\epsilon} \mid X\right) \\ &= (X'X)^{-1}X' \operatorname{Var}\left( \boldsymbol{\epsilon} \mid X\right) X(X'X)^{-1} \end{align*}

We assumed that that $\epsilon_i$ were IID. Hence for all $i$ we have $\operatorname{Var}(\epsilon_i) = \sigma^2$, and for the whole vector $\boldsymbol{\epsilon})$ we have $\operatorname{Var}(\boldsymbol{\epsilon}) = \sigma^2 I$ where $I$ is the identity matrix.

Continuing: \begin{align*} \operatorname{Var}\left( \hat{\mathbf{b}} \mid X\right) &= (X'X)^{-1}X' \sigma^2 IX(X'X)^{-1} \\ &= \sigma^2 (X'X)^{-1}X'X(X'X)^{-1} \\ &= \sigma^2 (X'X)^{-1} \end{align*}

$\sigma^2$ is a scalar and can be moved wherever by the commutative property of multiplication. $(X'X)^{-1}X'X = I$ by definition of an inverse matrix.

Let lowercase bold letters denote vectors. The linear model is:

$$ y_i = \mathbf{x}_i' \mathbf{b} + \epsilon_i $$

In matrix notation for all observations $i=1, \ldots, n$:

$$\mathbf{y} = X \mathbf{b} + \boldsymbol{\epsilon} $$

The OLS estimator for $\mathbf{b}$ is:

$$\hat{\mathbf{b}} = (X'X)^{-1}X'\mathbf{y} $$

Substituting:

$$\hat{\mathbf{b}} = (X'X)^{-1}X' \left( X \mathbf{b} + \boldsymbol{\epsilon}\right) $$

What's the variance of our estimator $\hat{\mathbf{b}}$? Take the variance (conditional on X)$:

\begin{align*} \operatorname{Var}\left( \hat{\mathbf{b}} \mid X\right) &= \operatorname{Var}\left( (X'X)^{-1}X' \left( X \mathbf{b} + \boldsymbol{\epsilon}\right) \mid X\right) \\ &= \operatorname{Var}\left( (X'X)^{-1}X' \boldsymbol{\epsilon} \mid X\right) \\ &= (X'X)^{-1}X' \operatorname{Var}\left( \boldsymbol{\epsilon} \mid X\right) X(X'X)^{-1} \end{align*}

We assumed that that $\epsilon_i$ were IID. Hence for all $i$ we have $\operatorname{Var}(\epsilon_i) = \sigma^2$, and for the whole vector $\boldsymbol{\epsilon})$ we have $\operatorname{Var}(\boldsymbol{\epsilon}) = \sigma^2 I$ where $I$ is the identity matrix.

Continuing: \begin{align*} \operatorname{Var}\left( \hat{\mathbf{b}} \mid X\right) &= (X'X)^{-1}X' \sigma^2 IX(X'X)^{-1} \\ &= \sigma^2 (X'X)^{-1}X'X(X'X)^{-1} \\ &= \sigma^2 (X'X)^{-1} \end{align*}

$\sigma^2$ is a scalar and can be moved wherever by the commutative property of multiplication. $(X'X)^{-1}X'X = I$ by definition of an inverse matrix.

You can estimate $\sigma^2$ using the residuals from your OLS regression. Let $e_i$ be the residual (as opposed to error-term) for observation $i$. That is:

$$ e_i = y_i - \mathbf{x}_i' \hat{\mathbf{b}} $$

Then the usual estimator for $\sigma^2$ is: $$\hat{\sigma}^2 = \frac{1}{n-k} \sum_i e_i^2$$ where $k$ is your number of regressors (including a constant).