#
# Written by:
# -- 
# John L. Weatherwax                2009-04-21
# 
# email: wax@alum.mit.edu
# 
# Please send comments and especially bug reports to the
# above email address.
# 
#-----

source('est_theta1_MA1.R')
source('est_phi1_theta1_AR1MA1.R')
source('est_theta1_theta2_MA2.R')

# Part (b) MA(1) model: 
#
r1 = -0.41
theta1 = est_theta1_MA1(r1)
print(sprintf("Part (b): theta1 = %10.6f", theta1))

# Part (c) ARIMA(1,0,1) model: 
#
r1 = 0.4
r2 = 0.32
phi_theta = est_phi1_theta1_AR1MA1(r1,r2) 
phi1 = phi_theta[1]
theta1 = phi_theta[2] # Root is theta_1 = 0.5
print(sprintf("Part (c): phi1 = %10.6f", phi1))
print(sprintf("Part (c): theta1 = %10.6f", theta1))


# Part (d) ARIMA(0,2,2) model: 
#
r1 = 0.62
r2 = 0.13
theta1_theta2 = est_theta1_theta2_MA2(r1,r2) # Root is theta_1 = -1.0803349, theta_2 = -0.2928770
theta1 = theta1_theta2[1]
theta2 = theta1_theta2[2] 
print(sprintf("Part (d): theta1 = %10.6f", theta1))
print(sprintf("Part (d): theta2 = %10.6f", theta2))


# Part (e) AR(2) model: 
#
r1 = 0.93
r2 = 0.81
phi1 = r1 * ( 1 - r2 )/( 1 - r1^2 ) 
phi2 = ( r2 - r1^2 )/( 1 - r1^2 ) 
print(sprintf("Part (e): phi1 = %10.6f", phi1))
print(sprintf("Part (e): phi2 = %10.6f", phi2))


# Do the calculations for Problem 6.2 here:
#
n = 101 
c0 = 0.25
se = sqrt(c0 * ( 1 + r1 ) * ( 1 - 2*r1^2 + r2 ) / ( n * (1-r1) * (1-r2) ) )
print(sprintf("Prob 6.2: hat{sigma}(bar{w})= %10.6f", se))

theta0 = ( 1 - 1.3079 + 0.40636 )* ( 0.23 ) 
print(sprintf("Prob 6.2: theta0= %10.6f", theta0))