graphics.off();
objects();
remove(list=objects());
options(scipen=20);
sink(file = "output.txt", append = FALSE, type = "output", split = TRUE);

library(mvtnorm);

DCC.LL <- function(theta, x)
{
lgl = 0;

# initialize vectors/matrices
alpha  <- 0;
beta   <- 0;
kappa  <- matrix(0, k, 1);
lambda <- matrix(0, k, 1);
Dsq    <- matrix(0, k, 1);
Q      <- matrix(0, k, k);
RR     <- matrix(0, k, k);
mu     <- matrix(0, 1, k);
epsL   <- matrix(0, k, 1);
rL     <- matrix(0, k, 1);
DsqL <- Dsq;
QL   <- Q;


# CHECK THAT PARAMETERS ADHERE TO RESTRICTIONS
i=1;
while(i<=npar){
  if(theta[i]<0){
    print("theta[i]<0 !");
    print(theta);
    return(Inf);
  }
  i=i+1;
}

# store variables into local variables kappa, lambda, alpha, and beta
# check for stationarity conditions
i=1;
while(i<=k){

  if(i==1){
    kappa[i] <- theta[i];
  }else{
    kappa[i] <- theta[(i*2)-1];
  }

  if(i==1){
    lambda[i] <- theta[i+1];
  }else{
    lambda[i] <- theta[(i*2)];
  }

  if(kappa[i]+lambda[i] >= 1){
    print("kappa[i] + lambda[i] >= 1 !");
    print(kappa[i]+lambda[i]);
    return(Inf);
  }
  i=i+1;
}

alpha <- theta[npar-1];
beta  <- theta[npar];

if(alpha+beta >= 1){
  print("alpha + beta >= 1 !");
  print(alpha+beta);
  return(Inf);
}

# -----------------------------------------------------------------------
# MAIN LOOP THAT CREATES LOGL

i=1;
while(i<=iend){

# create conditional variances
  if(i == 1){
    Dsq <- uncondv;
  }else{	
    rL <- x[i-1,1:k];
    j=1;
    while(j<=k){
      Dsq[j] <- uncondv[j]*(1-kappa[j]-lambda[j])+(kappa[j]*rL[j]^2)+(lambda[j]*DsqL[j]);
      j=j+1;
    }
  }	

# create conditional correlations
  if(i == 1){
    Q <- corr;
  }else{
    epsL <- t(x[i-1,1:k])/sqrt(DsqL);
    j=1;
    while(j<=k){
      p=j;
      while(p<=k){
        Q[p,j] <- corr[p,j]*(1-alpha-beta)+(alpha*epsL[p]*epsL[j])+beta*QL[p,j];
        Q[j,p] <- Q[p,j];
        p=p+1;
      }	
      j=j+1;
    }
  }	
  
# Form H matrix
  j=1;
  while(j<=k){
    p=j;
    while(p<=k){
      RR[p,j] <- Q[p,j]/sqrt(Q[p,p]*Q[j,j]);
      RR[j,p] <- RR[p,j];
      p=p+1;
    }
    j=j+1;
  }
  j=1;

  while(j<=k){
    p=j;
    while(p<=k){
      H[p,j,i] <<- RR[p,j]*sqrt(Dsq[p]*Dsq[j]);
      H[j,p,i] <<- H[p,j,i];
      p=p+1;
    }
    j=j+1;
  }
  
# Calculate LOGL

  lgl = lgl + dmvnorm(x[i,1:k], mu ,H[1:k,1:k,i] , log=TRUE);

  
#  A    = as.matrix(H[1:k,1:k,i]);
#  Ainv = solve(A);
#  detA = det(A);
#  sm   = (x[i,1:k]-mu)%*%Ainv%*%t(x[i,1:k]-mu);
#  logl = logl -(k*.5*log(2*pi)) -(.5*log(detA)) -(.5*sm);
  DsqL <- Dsq;
  QL   <- Q;
  
  i=i+1;
}

# Check that LOGL is in reasonable range and return value to optimizer
if(abs(lgl) >= 1.e-300 && abs(lgl) <= 2.e300){
                                        #      print(theta);
                                        #      print(logl);
                                        #      print(b);
  b <<- b+1;
  return(-lgl);
}else{
  print(".");
  print("log-likelihood outside allowed range!");
  print(".");
  return(Inf);


} }




# MAIN PROGRAM


myData <- read.table("data3.txt", header=FALSE);

ret    <- as.matrix(myData);
ret <- ret[,1:3];

ret    <- ret * 100;
k      <- ncol(ret);

nobs   <- nrow(ret);
print(k);
print(nobs);
iend   <- 900;
npar   <- 2*k + 2;
H      <- array(0, dim = c(k, k, iend));
   # var(ret[1:iend, 3]), var(ret[1:iend, 4]),var(ret[1:iend, 5]),var(ret[1:iend, 6])
uncondv <- diag( cov(ret) );
#  c(var(ret[1:iend, 1]), var(ret[1:iend, 2]));

corr    <- cor(ret[1:iend, 1:k]);
theta      <- c(0.09906,0.8967,0.047,0.94,0.04,0.94,0.006,0.98);

b <- 1;


DCC.LL(theta, x=ret);

#result   <- nlm(DCC.LL, theta, x=ret, hessian=T, steptol=.01, iterlim=10000,print.level=2);
result   <- nlm(DCC.LL, theta, x=ret, hessian=T, steptol=.00001, iterlim=10000,print.level=2);
#print(result);


invhess <- solve( result$hessian );
print( cbind(result$estimate,sqrt( diag(invhess) ),result$estimate/sqrt( diag(invhess) ) ) );


corrH      <- array(0, dim = c(k, k, iend));

ii <- 1;
while( ii <= iend){

  j=1;
  while(j<=k){
    p=j;
    while(p<=k){
      corrH[p,j,ii] <- H[p,j,ii]/sqrt(H[p,p,ii]*H[j,j,ii]);
      corrH[j,p,ii] <- corrH[p,j,ii];
      p=p+1;
    }
    j=j+1;
  }
#  print("H and corrH");
#  print( H[,,ii]);
#  print( corrH[,,ii]);
  ii=ii+1;
}


write.table(cbind(corrH[1,2,],corrH[1,3,],corrH[2,3,]) ,
            "corr.txt",sep=" ",row.names=FALSE,col.names = FALSE);




sink(file = NULL);