#
# 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.
# 
#-----

set.seed(0)

# Exercise 3 (EPage 411)
# 
zero_prices = read.table('../../BookCode/Chapter_14/ZeroPrices.txt',header=TRUE)
attach(zero_prices)

# Example code for fitting on EPage 384 / analytic form for the yield curve given on EPage 383
#

# To get good estimate of the coefficients we need to perhaps fit this model sequentially: 
#

# First: get an initially good estimate for theta0:
# 
fit_Nelson_Siegel_pt1 = nls(
  price ~ 100 * exp( - theta0 * maturity ),
  data=zero_prices,
  start=list(theta0=0.01),
  control=list(minFactor=1.e-5,warnOnly=TRUE) )
theta0_initial_estimate = summary(fit_Nelson_Siegel_pt1)$coef[,1] 

# Second: using the above estimate of theta0 get good estimate of theta2 and theta3 (and perhaps modify theta0 some):
# 
fit_Nelson_Siegel_pt2 = nls(
  price ~ 100 * exp( - ( theta0 + (theta2/theta3) * (1 - exp(-theta3*maturity)) / (theta3*maturity) - (theta2/theta3)*exp(-theta3*maturity) ) * maturity ),
  data=zero_prices,
  start=list(theta0=theta0_initial_estimate,theta2=0.5,theta3=0.5),
  control=list(minFactor=1.e-5,warnOnly=TRUE) )
coefs = summary(fit_Nelson_Siegel_pt2)$coef[,1] 
theta0_initial_estimate = coefs[1]
theta2_initial_estimate = coefs[2]
theta3_initial_estimate = coefs[3]

# Third: using the above estimates of theta0, theta2, and theta3 get an estimate of theta1 (perhaps modifying the other thetas some): 
#
fit_Nelson_Siegel = nls(
  price ~ 100 * exp( - ( theta0 + (theta1 + theta2/theta3) * (1 - exp(-theta3*maturity)) / (theta3*maturity) - (theta2/theta3)*exp(-theta3*maturity) ) * maturity ),
  data=zero_prices,
  start=list(theta0=theta0_initial_estimate,theta1=0.,theta2=theta2_initial_estimate,theta3=theta3_initial_estimate),
  control=list(minFactor=1.e-5,warnOnly=TRUE) )
coefs = summary(fit_Nelson_Siegel)$coef[,1]

# Plot our model with the data:
#
save_plots = F

if( save_plots ){ postscript("../../WriteUp/Graphics/Chapter14/chap_14_ex_3_yield_curve_fit.eps", onefile=FALSE, horizontal=FALSE) } 

plot( zero_prices$maturity, zero_prices$price )
theta0 = coefs[1]
theta1 = coefs[2]
theta2 = coefs[3]
theta3 = coefs[4] 
P_hat = 100 * exp( - ( theta0 + (theta1 + theta2/theta3) * (1 - exp(-theta3*maturity)) / (theta3*maturity) - (theta2/theta3)*exp(-theta3*maturity) ) * maturity )
lines( maturity, P_hat, type="l", col="red" )

if( save_plots ){ dev.off() }