##                                        
## maximum de vraisemblance avec Vasicek
##                                       
library("xtable")

## importation des données
data <- read.csv("usgg.csv",header=TRUE, sep=";")[(1:50)*30,]
delta <- 1/12
rf <- data$usgg3m/100
rf[which(rf<=0)] <- 0.001

##
## graphique quantile quantile
##
pdf("MLE-qqplot-norm.pdf")
qqnorm(diff(rf))
qqline
dev.off()

## fonction objectif de log-vraisemblance négatif

VASICEK.FUNOBJ <- function(PAR,RF,DELTA)
  {
    mu <- PAR[2]
    alpha <- PAR[1]
    sigma <- PAR[3]
    n <- length(RF)
    moyenne <- RF[-n] * exp(-alpha*delta)+mu*(1-exp(-alpha*delta))
    variance <- sigma^2 * (1-exp(-2*alpha*delta)) / (2*alpha)
    -sum(dnorm(RF[-1],moyenne,sqrt(variance),log=TRUE))
  }

ES1 <- optim(c(.5,-.5,.5),fn=VASICEK.FUNOBJ,RF=rf,DELTA=delta)

##                                        
## maximum de vraisemblance avec CIR (quasi-vraisemblance normale)
##

CIR.QV.FUNOBJ <- function(PAR,RF,DELTA)
  {
    mu <- PAR[2]
    alpha <- PAR[1]
    sigma <- PAR[3]
    n <- length(RF)
    moyenne <- RF[-n] * exp(-alpha * DELTA) +
      mu * (1-exp(-alpha * DELTA))
    variance <- RF[-n] * sigma^2/alpha * (exp(-alpha * DELTA) - exp(-2*alpha * DELTA)) +
      mu * sigma^2/2/alpha*(1-exp(-alpha * DELTA))^2
    -sum(dnorm(RF[-1],moyenne,sqrt(variance),log=TRUE))
  }

ES3 <- optim(c(alpha.st,mu.st,sigma.st),fn=CIR.QV.FUNOBJ,RF=rf,DELTA=delta)

##
## sorties
##
ES1$par[3] <- ES1$par[3]^2
ES3$par[3] <- ES3$par[3]^2

sink("MLE-dates.tex",append=FALSE,split=FALSE)
as.vector(data$date[c(1,length(data$date))])
sink()

sink("MLE-param.tex",append=FALSE,split=FALSE)
  xtable(cbind(Param=c("alpha","mu","sigma"),
               Vasicek=ES1$par,CIR.QL=ES3$par),
         caption="Paramètres estimés par maximum de vraisemblance pour 2 modèles",
         label="tab:estimParam")
sink()

ES1$par
ES3$par

##
## test de ratio de vraisemblance
##

loglik <- c(ES1$value,ES3$value)
n <- length(rf)-1

1-pchisq(-2*n*(loglik[2]-loglik[1]),4)