#
# Load the individual series and estimate parameters for each of the proposed models
#
# 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_theta1_theta2_MA2.R')
source('est_phi1_theta1_AR1MA1.R')

source('wt_std_error_AR1.R')
source('wt_std_error_MA1.R')
source('wt_std_error_MA2.R')
source('wt_std_error_AR1MA1.R')

# Series A: 
# 

Z = scan("../../Data/series_a.dat",strip.white=T)

# The ARIMA(1,0,1) model:
#
r = acf( Z, plot=F ); r1 = r[[1]][[2]]; r2 = r[[1]][[3]] 
phi1_theta1 = est_phi1_theta1_AR1MA1(r1,r2)
phi1 = phi1_theta1[1]
theta1 = phi1_theta1[2]

zbar = mean(Z) # do we need to add a mean?
c0 = var(Z) 
sigma_w_bar = wt_std_error_AR1MA1(length(Z),c0,r1,r2)

theta0 = ( 1. - phi1 ) * zbar
sigmaa2 = ( ( 1 - phi1^2 )/( 1 + theta1^2 - 2 * phi1 * theta1 ) ) * c0 

print(sprintf("Series A ARIMA(1,0,1): wbar= %10.6f; sigma_hat_wbar= %10.6f; c0= %10.6f", zbar,sigma_w_bar,c0))
print(sprintf("Series A ARIMA(1,0,1): phi1= %10.6f; theta1 = %10.6f; theta0= %10.6f; sigmaa2= %10.6f", phi1,theta1,theta0,sigmaa2))

# The ARIMA(0,1,1) model: 
#
Zd = diff(Z)
r = acf( Zd, plot=F ); r1 = r[[1]][[2]]; r2 = r[[1]][[3]]
theta1 = est_theta1_MA1(r1)

zbar = mean(Zd)
c0 = var(Zd)
sigma_w_bar = wt_std_error_MA1(length(Zd),c0,r1,r2)

theta0 = zbar
sigmaa2 = c0 / ( 1 + theta1^2 )

print(sprintf("Series A ARIMA(1,1,0): wbar= %10.6f; sigma_hat_wbar= %10.6f; c0= %10.6f", zbar,sigma_w_bar,c0))
print(sprintf("Series A ARIMA(1,1,0): theta1= %10.6f; theta0 (not significant)= %10.6f; sigmaa2= %10.6f", theta1,theta0,sigmaa2))


# Series B: 
# 

Z = scan("../../Data/series_b.dat",strip.white=T)

# The ARIMA(0,1,1) model: 
#
Zd = diff(Z)
r = acf( Zd, plot=F ); r1 = r[[1]][[2]]; r2 = r[[1]][[3]]
theta1 = est_theta1_MA1(r1)

zbar = mean(Zd)
c0 = var(Zd)
sigma_w_bar = wt_std_error_MA1(length(Zd),c0,r1,r2)

theta0 = zbar
sigmaa2 = c0 / ( 1 + theta1^2 )

print(sprintf("Series B ARIMA(0,1,1): wbar= %10.6f; sigma_hat_wbar= %10.6f; c0= %10.6f", zbar,sigma_w_bar,c0))
print(sprintf("Series B ARIMA(0,1,1): theta1= %10.6f; theta0 (not significant)= %10.6f; sigmaa2= %10.6f", theta1,theta0,sigmaa2))


# Series C: 
# 

Z = scan("../../Data/series_c.dat",strip.white=T)

# The ARIMA(1,1,0) model: 
#
Zd = diff(Z)
r = acf( Zd, plot=F ); r1 = r[[1]][[2]]; r2 = r[[1]][[3]]
phi1 = r1

zbar = mean(Zd)
c0 = var(Zd)
sigma_w_bar = wt_std_error_AR1(length(Zd),c0,r1,r2)

theta0 = zbar * ( 1 - phi1 )
sigmaa2 = c0 * ( 1 - phi1 * r1 )

print(sprintf("Series C ARIMA(1,1,0): wbar= %10.6f; sigma_hat_wbar= %10.6f; c0= %10.6f", zbar,sigma_w_bar,c0))
print(sprintf("Series C ARIMA(1,1,0): phi1= %10.6f; theta0 (not significant)= %10.6f; sigmaa2= %10.6f", phi1,theta0,sigmaa2))

# The ARIMA(0,2,2) model: 
#
Zd = diff(Z,differences=2)
r = acf( Zd, plot=F ); r1 = r[[1]][[2]]; r2 = r[[1]][[3]]
theta1_theta2 = est_theta1_theta2_MA2(r1,r2)
theta1 = theta1_theta2[1]
theta2 = theta1_theta2[2]

zbar = mean(Zd)
c0 = var(Zd)
sigma_w_bar = wt_std_error_MA2(length(Zd),c0,r1,r2)

theta0 = zbar
sigmaa2 = c0 / ( 1 + theta1^2 + theta2^2 )

print(sprintf("Series C ARIMA(0,2,2): wbar= %10.6f; sigma_hat_wbar= %10.6f; c0= %10.6f", zbar,sigma_w_bar,c0))
print(sprintf("Series C ARIMA(0,2,2): theta1= %10.6f; theta2= %10.6f; theta0 (not significant)= %10.6f; sigmaa2= %10.6f", theta1,theta2,theta0,sigmaa2))


# Series D: 
# 

Z = scan("../../Data/series_d.dat",strip.white=T)

# The ARIMA(1,0,0) model:
#
r = acf( Z, plot=F ); r1 = r[[1]][[2]]; r2 = r[[1]][[3]] 
phi1 = r1

zbar = mean(Z) # do we need to add a mean?
c0 = var(Z) 
sigma_w_bar = wt_std_error_AR1(length(Z),c0,r1,r2)

theta0 = ( 1. - phi1 ) * zbar
sigmaa2 = ( 1 - phi1*r1 ) * c0 

print(sprintf("Series D ARIMA(1,0,0): wbar= %10.6f; sigma_hat_wbar= %10.6f; c0= %10.6f", zbar,sigma_w_bar,c0))
print(sprintf("Series D ARIMA(1,0,0): phi1= %10.6f; theta0= %10.6f; sigmaa2= %10.6f", phi1,theta0,sigmaa2))

# The ARIMA(0,1,1) model: 
#
Zd = diff(Z)
r = acf( Zd, plot=F ); r1 = r[[1]][[2]]; r2 = r[[1]][[3]]
theta1 = est_theta1_MA1(r1)

zbar = mean(Zd)
c0 = var(Zd)
sigma_w_bar = wt_std_error_MA1(length(Zd),c0,r1,r2)

theta0 = zbar
sigmaa2 = c0 / ( 1 + theta1^2 )

print(sprintf("Series D ARIMA(0,1,1): wbar= %10.6f; sigma_hat_wbar= %10.6f; c0= %10.6f", zbar,sigma_w_bar,c0))
print(sprintf("Series D ARIMA(0,1,1): theta1= %10.6f; theta0 (not significant)= %10.6f; sigmaa2= %10.6f", theta1,theta0,sigmaa2))



## Z = scan("../../Data/series_e.dat",strip.white=T)

## Z = scan("../../Data/series_f.dat",strip.white=T)