高斯一牛顿迭代法(Gauss-Newton iteration method) 是非线性回归模型中求回归参数进行最小二乘的一种迭代方法，该法使用泰勒级数展开式去近似地代替非线性回归模型，然后通过多次迭代，多次修正回归系数，使回归系数不断逼近非线性回归模型的最佳回归系数，最后使原模型的残差平方和达到最小。 本文使用Gauss-Newton算法去求解时间序列模型的参数，一个ARMA(1,1)，一个AR(1).

1.

使用gauss-newton algorithm去拟合ARMA(1,1). 找出\(\phi\), \(\theta\), \(\sigma_w^2\)这三个参数.

library(tidyverse)
## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --
## v ggplot2 3.3.5     v purrr   0.3.4
## v tibble  3.1.4     v dplyr   1.0.7
## v tidyr   1.1.3     v stringr 1.4.0
## v readr   2.0.1     v forcats 0.5.1
## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::first()  masks xts::first()
## x dplyr::lag()    masks stats::lag()
## x dplyr::last()   masks xts::last()
dat=read.table("https://canvas.illinois.edu/files/2382460/download?download_frd=1&verifier=Q3JkKvukvxw2IkPSg9J5PShPbJSArRm95RKqqibq")
head(dat)
##            x
## 1 -4.1845303
## 2 -4.1398843
## 3 -3.6553543
## 4 -2.1845921
## 5 -0.5192022
## 6 -0.9634690

\[X_t-\phi X_{t-1}=W_t+\theta W_{t-1}\]

\[\rho(h)=\frac{(1+\theta \phi)(\phi+\theta)}{\phi(1+2\theta\phi+\theta^2)}\phi^h\]

x = dat$x # data
r = acf1(x, 1, plot=FALSE) # acf(1)

num = length(x) 
## Estimation

niter = 10000
para=matrix(rep(0,2*(niter+1)),nrow = niter+1, ncol = 2)
para-> Sz -> Szw
sc=c(0)
c(0)->w 

matrix(rep(0,2*num),nrow = num,ncol = 2)-> z
para[1,] = c(1,1) # MME
for (j in 1:niter){
  Szz=matrix(c(0,0,0,0),nrow = 2)
  for (i in 2:num){ 
  w[i] = x[i] - para[j,1]*x[i-1]-para[j,2]*w[i-1]
  z[i,2] = w[i-1] - para[j,2]*z[i-1,2]
  z[i,1] = x[i-1] - para[j,2]*z[i-1,1]
  Szz=Szz+z[i,]%*%t(z[i,])
  }
  sc[j]=sum(w*w)
  Szw[j,] = t(colMeans(z*w))
  para[j+1,] = as.double(para[j,] + t(solve(Szz)%*%Szw[j,]))
}
# Results
output=cbind(iteration=1:niter-1, phi=para[1:niter,1], thetahat=para[1:niter,2], sc=sc[1:niter])
plot(output[,2])

print(output[niter,2])
##      phi 
## 0.939893
plot(output[,3])

print(output[niter,3])
##  thetahat 
## 0.5412998

As we can see above, the estimation of \(\phi\) is 0.939893, \(\theta\) is 0.5412998.

When it comes to \(\sigma^2\)

print(output[niter,4])
##       sc 
## 184.2437

\[\sigma^2=\frac{184.2437}{200-1}=0.9258477\]

Compare with the function provided by R.

arima(dat$x, order=c(1,0,1), method="CSS")
## 
## Call:
## arima(x = dat$x, order = c(1, 0, 1), method = "CSS")
## 
## Coefficients:
##          ar1     ma1  intercept
##       0.9019  0.5511    -2.6397
## s.e.  0.0315  0.0604     1.0641
## 
## sigma^2 estimated as 0.9092:  part log likelihood = -274.27

While the estimation of \(\phi\) is 0.9019, \(\theta\) is 0.5511 and \(\sigma^2\) is 0.9092. It’s approximately the same.

2.

理论推导出AR(1)的gauss-newton algorithm的估计。

\[X_t = \phi X_{t-1}+W_t\]

We can write the error as

\[w_t(\phi)=x_t-\phi x_{t-1}(\phi)\]

So our goal is to find a value of \(\phi\) to minimize \[S_c(\phi)=\sum_{t=1}^n w_t^2(\phi)\] let \(\phi_{(0)}\) be the initial estimate of \(\phi\). And then we use a first-order Taylor approximation of \(w_t\) at \(\phi_{(0)}\).

\[S_c(\phi)=\sum_{t=1}^n w_t^2(\phi)\approx \sum_{t=1}^n[w_t(\phi_{(0)})-(\phi-\phi_{(0)})z_t(\phi_{(0)})]^2\] Where

\[z_t(\phi_{(0)})=-\frac{\partial w_t(\phi)}{\partial \phi}|_{\phi=\phi_{(0)}}\]

as we can also find that

\[\frac{\partial w_t(\phi)}{\partial \phi}=-x_{t-1}(\phi)\]

We set \(\frac{\partial w_0(\phi)}{\partial \phi}=0\)

So \[z_t(\phi)=x_{t-1}(\phi)\] \[Q(\phi)=\sum_{t=1}^n[w_t(\phi_{(0)})-(\phi-\phi_{(0)})x_{t-1}(\phi_{(0)})]^2\] We should minimize this function

\[\hat{\phi}-\phi_{(0)}=\frac{\sum_{t=1}^n x_{t-1}(\phi_{(0)})w_t(\phi_{(0)})}{\sum_{t=1}^nx_{t-1}^2(\phi_{(0)})}\]

Then we have

\[\hat{\phi}=\phi_{(0)}+\frac{\sum_{t=1}^n x_{t-1}(\phi_{(0)})w_t(\phi_{(0)})}{\sum_{t=1}^nx_{t-1}^2(\phi_{(0)})}\]

Finally, the Gauss-Newton procedure in this case is

\[\phi_{(j+1)}=\phi_{(j)}+\frac{\sum_{t=1}^n x_{t-1}(\phi_{(j)})w_t(\phi_{(j)})}{\sum_{t=1}^nx_{t-1}^2(\phi_{(j)})}\]

\[\phi_{(j+1)}=\phi_{(j)}+\frac{\sum_{t=1}^n x_{t-1}(\phi_{(j)})(x_t-\phi_{(j)} x_{t-1}(\phi_{(j)}))}{\sum_{t=1}^nx_{t-1}^2(\phi_{(j)})}\\=\frac{\sum_{t=1}^n x_{t-1}(\phi_{(j)})x_t}{\sum_{t=1}^nx_{t-1}^2(\phi_{(j)})}\]

It’s a conditional estimator.