Maximum Likelihood with dependent observations In the context of a normally distributed dependent variable, the likelihood for a single observation is given as:
$L(y_i|x,β,σ2)=\displaystyle\cfrac{1}{√2πσ^2}e^{\displaystyle\cfrac{(y_i−x_iβ)}{2σ^2}}$
Then, given that all observations are independent, the likelihood for the whole sample is: 
$L(y_i|x,β,σ2)= \displaystyle \prod_{i=1}^{\infty} \displaystyle\cfrac{1}{√2πσ^2}e^{\displaystyle\cfrac{(y_i−x_iβ)}{2σ^2}}$
In the matrix form, the log likelihood would be: 
$lnL = - {\displaystyle\cfrac{N}{2}}ln(2\pi) - {\displaystyle\cfrac{N}{2}}ln(σ^2)- {\displaystyle\cfrac{1}{2σ^2}}(y - X\beta)'(y - X\beta)$
My question is: how would you adjust both the likelihood and the log-likelihood formula to take into account the fact that observations are not independent? If you could explain the logic of how the formula is changed I would be immensely grateful! 
 A: As Whuber points out it depends on the assumption. So here is an example of some assumptions that could be made. First stack the individual equations
$$y_i = x_i^\top \beta + u_i$$
to get multivariate version
$$y = X\beta + u$$
then make some assumptions
$$\mathbb E[u] = 0$$
$$\mathbb Var(u) = \mathbb E[uu^\top] = \Omega $$
$$ det(\Omega) \not = 0$$
$$ rank(X) = K \ and \ X \ is\ fixed $$
$$u \sim Multivariate\ Normal$$
then the log likelihood can be written
$$\frac{1}{2} \log \det(\Omega^{-1}) - \frac{1}{2} (y - X\beta)^\top\Omega^{-1}(y - X\beta)$$
Assume that you want to estimate $\beta$ and that $\Omega$ is known, then maximizing the likelihood is the same as minimizing the sum of residuals weighted by the inverse of the variance
$$S(\beta) = (y - X\beta)^\top\Omega^{-1}(y - X\beta)$$
Of course, this gives rise to the Generalized Least Squares estimator, see GLS-Wiki, and the estimate is
$$\hat \beta_{GLS} = (X^\top \Omega^{-1} X)^{-1}X^\top \Omega^{-1}y.$$
$\Omega^{-1}$ is rarely known, but the GLS estimator is still theoretically important due to its efficiency properties.
If $\hat \Omega$ consistently estimates the covariance matrix, then the feasible GLS estimator is defined as
$$\hat \beta_{FGLS} = (X^\top \hat\Omega^{-1} X)^{-1}X^\top \hat \Omega^{-1}y$$
and if $\hat \beta_{OLS} := (X^\top X)^{-1}X^\top y,$ it is possible to get OLS residuals
$$\hat u = y - X \hat \beta_{OLS},$$ and then estimate $\Omega$ using these residuals. However, each term in the covariance matrix $\sigma_{ij}$ is only realized once in the residuals $\hat u_i \hat u_j,$ or, more importantly, this is not a consistent estimate of the variance, hence some structure must be imposed.
Amemiya (1985) Advanced Econometrics, page 185, discusses several standard approaches.
Other important sources imposing minimal assumptions on the covariance matrix are
White (1980) "A Heteroskedasticity-Consistent Covariance Matrix Estimator and a Direct Test for Heteroskedasticity" in Econometrica,
and
Newey, Whitney K; West, Kenneth D (1987). "A Simple, Positive Semi-definite, Heteroskedasticity and Autocorrelation Consistent Covariance Matrix". Econometrica.
For an explanation of implementation of these covariance estimators in R see
Zeileis, A. (2004). "Econometric Computing with HC and HAC Covariance Matrix Estimators". Journal of Statistical Software.
$$$$
