## com11c of march 2022,
#  MA(3), estimation of parameters 
#  x_t = w_t + \theta_1 w_{t-1} + \theta_2 w_{t-2} + \theta_3 w_{t-3} 

qq <- 3
theta0 <- 1
theta1 <- 0.55
theta2 <- 0.33
theta3 <- 0.11

sigmaw <- 1.00 

nn <- 10^4
xx <- 0*(1:nn)
ww <- rnorm(nn)

for (j in 4:nn)
{
xx[j] <- theta0*ww[j] + theta1*ww[j-1] + theta2*ww[j-2] + theta3*ww[j-3]
}

plot(xx,type="o")

rhohat1 <- acf(xx)$acf[2]
rhohat2 <- acf(xx)$acf[3]
rhohat3 <- acf(xx)$acf[4]

## then estimating rho1 and rho2 
#  from equation formulae with empirical values

Q3 <- function(para)
{
theta1 <- para[1]
theta2 <- para[2]
theta3 <- para[3]
# 
varr3 <- 1 + theta1^2 + theta2^2 + theta3^2 
rho1 <- (theta0*theta1 + theta1*theta2 + theta2*theta3)/varr3
rho2 <- (theta0*theta2 + theta1*theta3)/varr3
rho3 <- (theta0*theta3)/varr3
# 
well1 <- rhohat1 - rho1
well2 <- rhohat2 - rho2
well3 <- rhohat3 - rho3
well1^2 + well2^2 + well3^2 
}

Q3(c(theta1,theta2,theta3))

nils3 <- nlm(Q3,c(theta1,theta2,theta3),hessian=T)

thetahat <- nils3$estimate
sigmawhat <- sqrt( var(xx) / (1+thetahat[1]^2+thetahat[2]^2+thetahat[3]^2) )

check <- c(thetahat[1],thetahat[2],thetahat[3],sigmawhat)

print(round(check,4))