## 3 Mathematical Model
### Calculating the \(p\)-Order Gauss-Markov Model Parameters
#### 3.3.2 Calculating the Parameters  
After estimating the best order \(p\), we can estimate the parameters \(A_{i},i=1,..p\), and once \(A_{i}\) is estimated, \(\sigma_{\emptyset}^{2}\) can be estimated in turn. There are several methods to estimate the coefficients of the \(p\) -order Gauss-Markov model \(AR(p)\), such as the Markov Chain Monte Carlo, the least square optimization method, and the method of moments using Yule Walker equations. We chose to use the method of moments through the Yule Walker equations, due to its simplicity and the fact that it produces consistent estimators. Following [26], we provide the estimation of \(A_{isp}\) and \(\sigma_{\emptyset}\) as follows:  
Referring to Equation (3), we proceed to solve for the parameters of \(AR(p)\). Multiplying Equation (3) by \(X_{n-p}\), we get:  
\[X_{n-p}X_{n}\ =\sum_{i=1}^{p}A_{i}X_{n-p}X_{n-i}\ +\ X_{n-p}\partial_{n} \tag{5}\]  
Taking the expectation of Equation (5), we get:  
\[E(X_{n-p}X_{n})\ =E\left(\sum_{l=1}^{p}A_{i}X_{n-p}X_{n-i}\ +\ X_{n-p} \partial_{n}\right) \tag{6}\]  
\[E(X_{n-p}X_{n})\ =\sum_{i=1}^{p}A_{i}E\big{(}X_{n-p}X_{n-l}\big{)}\ \ \ +\ E(X_{n-p}\partial_{n}) \tag{7}\]  
\[E(X_{n-p}\partial_{n})=\begin{cases}0,p\
eq 0\\ \sigma_{\emptyset}^{2},p=0\end{cases} \tag{8}\]  
Equation (8) follows from the observation that \(X_{n-p}\) and \(\partial_{n}\) are uncorrelated for \(p\
eq 0\), because \(\partial_{n}\) of the current time is unrelated to \(X_{n-p}\) and thus uncorrelated with previous values of the process, hence,  
\[E(X_{n-p}X_{n})\ =\sum_{i=1}^{p}A_{i}E\big{(}X_{n-p}X_{n-i}\big{)} \tag{9}\]  
Dividing Equation (9) by N, the observation window of the process \(X_{n}\), and using the fact that the process is stationary, E(\(X_{n-p}X_{n}\)), is a function of the time lag, \(p\), only, and E(\(X_{n-p}X_{n}\)), can be denoted by \(c_{p}\), also using \(c_{p}=c_{-p}\), Equation (9) can be written as  
\[c_{p}=\ \sum_{l=1}^{p}A_{i}c_{i-p} \tag{10}\]  
Dividing \(c_{p}\) by \(c_{0}\), we get:\[r_{p}=\sum_{i=1}^{p}A_{i}r_{i-p} \tag{11}\]  
where \(r_{p}\) is the autocorrelation coefficient. Rewriting the equations for all \(r_{1\leq i\leq p}\); _i.e._, \(r_{1}=\sum_{i=1}^{p}A_{i}r_{i-1}\), \(r_{2}=\sum_{i=1}^{p}A_{i}r_{i-2}\), \(r_{3}=\sum_{i=1}^{p}A_{i}r_{i-3}\)...., \(r_{p}=\sum_{i=1}^{p}A_{i}r_{i-p}\) in matrix form, we have:  
\[R=\beta A \tag{12}\] \[R=\left[\begin{array}{c}r_{1}\\...\\ r_{p}\end{array}\right]\text{, }A=\left[\begin{array}{c}A_{1}\\...\\ A_{p}\end{array}\right]\text{, and }\beta=\left[\begin{array}{ccc}r_{0}&r_{1}&...&r_{p-1}\\ r_{1}&r_{0}&r_{p-2}\\ &\vdots&\vdots\\ r_{p-1}&r_{1}&r_{0}\end{array}\right] \tag{13}\]  
Note that \(\beta\) is full-rank and symmetric matrix, hence it is guaranteed that its inverse exist. Thus, \(A\) can be estimate from Equation (12) as \(\widetilde{A}=\beta^{-1}R\).  
Once the model parameters are estimated, the noise variance \(\sigma_{\emptyset}^{2}\) can be estimated by setting \(p=0\) in Equation (7):  
\[E(X_{n}X_{n})=\sum_{i=1}^{p}A_{i}E(X_{n}X_{n-i})\quad+\ E(X_{n}\partial_{n}) \tag{14}\]  
Using the fact that \(X_{n}\) is a stationary process again, \(c_{i}=\ c_{-i}\), and \(\sigma_{\emptyset}^{2}=E(X_{n}\partial_{n})\), Equation (14) can be written as:  
\[c_{0}=\sum_{i=1}^{p}A_{i}c_{i}+\sigma_{\emptyset}^{2} \tag{15}\] \[\sigma_{\emptyset}^{2}=c_{0}-\sum_{i=1}^{p}A_{i}c_{i} \tag{16}\]