###############################################################################################
################################UNIVARIATE CASE ################################################
###############################################################################################

DISTRnovasARCH_function(cc, a0 = 0.1, scale = 1)
{
# this is the distr of Xt in general NOVAS 
# i.e. Xt=scale*W/srqt(1-a0*W^2) 
# where W is N(0,1) TRUNCATED to range +- 1/sqrt(a0)
# and  scale=sqrt(alpha*varX+ sum_{i=1}^p a_i X_{t-i} )
# cc can vector of values
# when a0 tends 0, this distr tends to N(0,1)
# it looks like this distr has almost finite 2nd moment, but infinite 3rd and 4th
	cc <- cc/scale
	distr <- exp(( - cc^2)/(2 * (1 + a0 * cc^2)))/((1 + a0 * cc^2)^(3/2))
	distr <- distr/(scale * sqrt(2 * 3.1416) * (pnorm(1/sqrt(a0)) - pnorm(-1/sqrt(a0))))
	return(distr)
}
###############################################################################################

RDISTRnovasARCH_function(n, a0 = 0.1)
{
# generates n random variables from distribution f(.;a0,1)
	n3 <- floor(1.3 * n)  # corrected 
	Z <- rnorm(n3)
	c <- 1/sqrt(a0)
	W <- Z[Mod(Z) < c]
	U <- W/sqrt(1 - a0 * W^2)
	U <- U[1:n]
	U
}
###############################################################################################

MLEgarch_function(a0 = 0.1, A = 0.05, B = 0.9, C = 3e-006, xseries)
{
# a0,A,B,C are starting values for the numerical MLE
# xseries is data series
ee_1e-08
bbb <- nlminb(start = c(a0, A, B, C), obj = NegLogLikGARCH, lo =c(ee, ee, ee, ee),  
up = c(0.2, 1, 1, 1), x =xseries)
	return(bbb$parameters, bbb$message)	
# return the final a0,A,B,C values plus a message indicating whether convergence was reached
}

###############################################################################################
NegLogLikGARCH_function(Avector, x)
{
# negative of log (pseudo) likelihood of returns series x
# to be minimized by  nlminb  to find: Avector 
# this function is similar to NegLogLik which is for ARCH
# Now: Avector is the vector of four GARCH parameters
# Avector =c(a0,A,B,C) in the notation of NoVaSDISTR paper
# all a0,A,B,C are positive and A+B<1
	n <- length(x)
	a0 <- Avector[1]
	A <- Avector[2]
	B <- Avector[3]
	C <- Avector[4]
	p <- 50
	a <- c(1:p) * 0
	Aconst <- C/(1 - B)
	for(i in 1:p) {
		a[i] <- A * (B^(i - 1))
	}
	yy <- as.ts(x)
	kk <- 0 * yy^2 + Aconst
	for(i in 1:p) {
		kk <- kk + a[i] * lag(yy^2,  - i)
	}
#kk is basically the Sigma; length(kk)=n-p
	kk <- as.vector(kk)
	xs <- yy[(p + 1):n]/sqrt(kk)
	LL <- 0
	for(i in 1:(n - p)) {
		ax2 <- 1 + a0 * (xs[i])^2
		LL <- LL + (xs[i])^2/(2 * ax2) + (3/2) * log(ax2) + log(pnorm(1/sqrt(a0)) - pnorm(-1/sqrt(a0)))
		LL <- LL + log(sqrt(kk[i])) + log(sqrt(6.28))
	}
	LL <- LL/(n - p)
	LL
}
 
###############################################################################################
################################BIVARIATE CASE ################################################
###############################################################################################
 
DISTR2novasARCH_function(uu, a0 = c(0.1, 0.1), scale = c(1, 1), Rho = 0)
{
# returns the value of the bivariate density f_rho (cc;a0,scale) at bivariate point uu
# uu can NOT NOW be a vector of points
	uu <- uu/scale
	Scale <- scale[1] * scale[2]
	gu <- uu * 0
	guoveru <- gu
	C0 <- 1/sqrt(a0)	
	# normfactor below is the prob. that the normal r.v. is within \pm C0
	normfactor <- pmvnorm(c(C0[1], C0[2]), rho = Rho) - pmvnorm(c( - C0[1], C0[2]), rho
		 = Rho) - pmvnorm(c(C0[1],  - C0[2]), rho = Rho) + pmvnorm(c( - C0[1],  - C0[2]), rho = Rho)
	for(i in 1:2) {
		gu[i] <- uu[i]/sqrt(1 + a0[i] * uu[i]^2)
		guoveru[i] <- 1/sqrt(1 + a0[i] * uu[i]^2)
	}
	frho <- (dmvnorm(gu, rho = Rho) * (guoveru[1])^3 * (guoveru[2])^3)/(normfactor * Scale)
	frho <- max(frho, 0.000001) # because need the log of that later
return(frho)
}

###############################################################################################

RDISTR2novasARCH_function(n, a0 = c(0.1, 0.1), Rho = 0)
{
# generates n random variables (pairs) from bivariate density f_rho(.;a0,1)
# Note: if a0 is too big, the result may contain fewer than n variables
	n3 <- floor(1.3 * n)
	Z <- rmvnorm(n3, rho = Rho)
	cc <- c(1:n3)
	C0 <- 1/sqrt(a0)
	cc1 <- cc[Mod(Z[, 1]) < C0[1]]
	cc2 <- cc[Mod(Z[, 2]) < C0[2]]
	cc <- intersect(cc1, cc2)
	W <- Z[cc,  ]
	U <- W
	U[1] <- W[1]/sqrt(1 - a0[1] * W[1]^2)
	U[2] <- W[2]/sqrt(1 - a0[2] * W[2]^2)
	U <- U[1:n,  ]
	U
}


###############################################################################################

NegLogLik2GARCH_function(Avector, x)
{
# negative of log (pseudo) likelihood of BIVARIATE returns series x, i.e., x is n by 2 matrix
# this function is to be minimized by  nlminb  to find: Avector 
# Avector is the 9-variate vector of four + four GARCH parameters plus the correlation Rho 
	n <- length(x[, 1])
	a0 <- c(1:2)
	A <- c(1:2)
	B <- c(1:2)
	C <- c(1:2)
	a0[1] <- Avector[1]
	A[1] <- Avector[2]
	B[1] <- Avector[3]
	C[1] <- Avector[4]
	a0[2] <- Avector[5]
	A[2] <- Avector[6]
	B[2] <- Avector[7]
	C[2] <- Avector[8]
	rho <- Avector[9]
	p <- 50	# provided sample size n > > 50
      # approximate GARCH(1,1) by ARCH(p) with p=50 
	a <- c(1:p) * 0
	Aconst <- C[1]/(1 - B[1])
	for(i in 1:p) {
		a[i] <- A[1] * (B[1]^(i - 1))
	}
	yy1 <- as.ts(x[, 1])
	kk1 <- 0 * yy1^2 + Aconst
	for(i in 1:p) {
		kk1 <- kk1 + a[i] * lag(yy1^2,  - i)
	}
	kk1 <- as.vector(kk1)
	a <- c(1:p) * 0
	Aconst <- C[2]/(1 - B[2])
	for(i in 1:p) {
		a[i] <- A[2] * (B[2]^(i - 1))
	}
	yy2 <- as.ts(x[, 2])
	kk2 <- 0 * yy2^2 + Aconst
	for(i in 1:p) {
		kk2 <- kk2 + a[i] * lag(yy2^2,  - i)
	}
	kk2 <- as.vector(kk2)
	LL <- 0
	for(i in 1:(n - p)) {
	LL <- LL - log(DISTR2novasARCH(x[(i+p),], c(Avector[1], Avector[5]), scale = c(kk1[i], kk2[i]), Rho = Avector[9]))
	}
	LL
}
###############################################################################################
NegLogLik2.1GARCH_function(Avec, xx){
# to be optimized by nlminb over Avec= the first 4 parameters of Avector
# xx is n by 3 matrix whose first two columns consist of the BIVARIATE returns series x which is n by 2 matrix
# and the third column has as its first 5 entries the remaining 5 parameters of Avector
X_xx[,1:2]
AVECTOR_c(1:9)
AVECTOR[1:4]_Avec
AVECTOR[5:9]_xx[1:5,3]
LL_NegLogLik2GARCH(AVECTOR,X)
LL}

###############################################################################################
NegLogLik2.2GARCH_function(Avec, xx){
# to be optimized by nlminb over Avec= the second 4 parameters of Avector
# xx is n by 3 matrix whose first two columns consist of the BIVARIATE returns series x which is n by 2 matrix
# and the third column has as its first 5 entries the remaining 5 parameters of Avector (the first 4 and the last)
X_xx[,1:2]
AVECTOR_c(1:9)
AVECTOR[5:8]_Avec
AVECTOR[1:4]_xx[1:4,3]
AVECTOR[9]_xx[5,3]
LL_NegLogLik2GARCH(AVECTOR,X)
LL}

###############################################################################################
NegLogLik2.3GARCH_function(Avec, xx){
# to be optimized by nlminb over Avec= the last parameter (correlation) of Avector
# xx is n by 3 matrix whose first two columns consist of the BIVARIATE returns series x which is n by 2 matrix
# and the third column has as its first 8 entries the remaining 8 parameters of Avector (the first 8)
X_xx[,1:2]
AVECTOR_c(1:9)
AVECTOR[1:8]_xx[1:8,3]
AVECTOR[9]_Avec
LL_NegLogLik2GARCH(AVECTOR,X)
LL}


################################################################################################
################################################################################################

DIVCONMLE2garch_function(avector, xseries, LAMBDA=0.2,loops=20){
# divide-and-conquer profile numerical MLE in the bivariate time series case
# avector consists of starting values for the bivariate numerical MLE
# avector =c( a0,A,B,C for first series, a0,A,B,C for second series, and rho)
# xseries  is bivariate data series i.e., it is n by 2 matrix
# NEEDS ACCURATE STARTING VALUES FOR THE FIRST 8 COORDINATES OF avector (e.g. obtained by univariate MLE)
# ALSO: SHOULD BE RUN WITH DIFFERENT STARTING VALUES AND THE LogLikelihood ENTRIES COMPARED
# convergence factor LAMBDA is in (0,1)
# loops is number of iterations; should be an integer in [3,100]
        n <- length(xseries[, 1])
        AVECTOR <- matrix(0, 9, 100)
        LogLikelihood_c(1:100)*0
        LogLikelihood[1]_ -NegLogLik2GARCH(avector, x = xseries)        
        LogLikelihood[2]_ LogLikelihood[1]
        AVECTOR[, 1] <- c(1:9) * 0 + 1
        AVECTOR[, 2] <- avector
#       anorm <- c(1:100) * 0
#       anorm[1] <- 9
#       anorm[2] <- sum((AVECTOR[, 2])^2)
 xxx <- cbind(xseries, c(1:n) * 0)
                        xxx[1:8, 3] <- avector[1:8]
aka_acf(xseries, lag.max=1, plot=F)
aka_aka[[1]]
   avector[9] <- aka[1,2,1] 
 # get starting value for the correlation rho as cross-correlation of series
#bb3 <- nlminb(start = avector[9], obj = NegLogLik2.3GARCH, lo = -0.99999 , up = 0.99999 , xx = xxx)
      #                  bb3 <- as.vector(bb3$parameters)
       #                 avector[9] <- bb3[1]    # get starting value for the correlation rho
                        #
                        for(k in 3:loops) {  # iteration starts 
avec_AVECTOR[,(k-1)]
lowerlim <- avec  * (1-.3*(1-k/101)) # get narrower limits as iteration proceeds
upperlim <- avec  * (1+.5*(1-k/101))
            for(i in 1:9) {     if(upperlim[i] >= 1) {  upperlim[i] <- 0.99999          }       }
        if(lowerlim[9] <=  -1) {        lowerlim[9]  _  -0.99999 }
                AVECTOR[, k] <-  AVECTOR[,(k-1)] 
                aaa_max(0.0001, sum ( ( AVECTOR[,(k-2)])^2))
       aaa_sum ( (AVECTOR[,(k-1)] -AVECTOR[,(k-2)])^2) / aaa
                if(aaa > 1e-006) { # loop as long as convergence is NOT reached 
                      #  avec_AVECTOR[,(k-1)]
                        xxx <- cbind(xseries, c(1:n) * 0)
                        xxx[1:4, 3] <- avec[1:4]
                        xxx[5, 3] <- avec[9]
                        bb2 <- nlminb(start = avector[5:8], obj = NegLogLik2.2GARCH,
gradient=NULL, hessian=NULL, scale=1, control=nlminb.control( iter.max = 100, rel.tol = .00001), 
lo = lowerlim[5:8], up = upperlim[5:8], xx = xxx)
                        bb2 <- as.vector(bb2$parameters)
                        avector[5:8] <- bb2      
                        # 
                        xxx <- cbind(xseries, c(1:n) * 0)
                        xxx[1:5, 3] <- avec[5:9]
bb1 <- nlminb(start = avector[1:4], obj = NegLogLik2.1GARCH, 
gradient=NULL, hessian=NULL, scale=1, control=nlminb.control( iter.max = 100, rel.tol = .00001), 
lo = lowerlim[1:4], up = upperlim[1:4], xx = xxx)
                        bb1 <- as.vector(bb1$parameters)
                        avector[1:4] <- bb1     
                        #
                        xxx <- cbind(xseries, c(1:n) * 0)
                        xxx[1:8, 3] <- avec[1:8]
bb3 <- nlminb(start = avector[9], obj = NegLogLik2.3GARCH, lo = lowerlim[9], up = upperlim[9], xx = 
xxx)
                        bb3 <- as.vector(bb3$parameters)
                        avector[9] <- bb3[1]            
                        AVECTOR[, k] <- LAMBDA*avector +(1-LAMBDA)*AVECTOR[,(k-1)]
                     avector_ AVECTOR[, k]      }
                    LogLikelihood[k]_ -NegLogLik2GARCH(avector, x = xseries)    }
           AVECTOR[, 1] _ AVECTOR[, 2] 
        rbind(AVECTOR,LogLikelihood) } 


###############################################################################################
###############################################################################################
##### HOW TO INTERACTIVELY USE THE DIVCONMLE2garch TO GET THE MLEs IN THE BIVARIATE CASE #####
# the data at hand is a bivariate time series: 1st coordinate series is series1, and 2nd is series2
aa_c(1:9)*0  # this will serve as starting values for DIVCONMLE2garch
bb_MLEgarch(xseries = series1)
bb$message # if bb$message indicates relative convergence was achieved proceed; otherwise, repeat ab
ove with careful starting values
aa[1:4]_bb$parameters
#
bb_MLEgarch(xseries = series2)
bb$message # if bb$message indicates relative convergence was achieved proceed; otherwise, repeat ab
ove with careful starting values
aa[5:8]_bb$parameters
#
par(mfrow=c(5,2)) # for plotting 
cc_DIVCONMLE2garch(aa,cbind(series1,series2)) # this step may take a while...
for (i in 1:10){ tsplot(cc[i,2:loops])}               # view plots to gauge convergence
# first 9 plots: evolution of the 9 parameter estimates over iterations--hopefully, showing converge
nce.
# The 9 parameters are: a0,A,B,C for first series, a0,A,B,C for second series, and rho.
# 10th plot: evolution of the value of the LogLikelihood correspodning to above 9 parameters--ideall
y, this is a
# (roughly) increasing graph that tapers off (converges) towards the right end showing the maximizat
ion is working.
# If all works as planned, the parameter MLEs are the final values, i.e.,  cc[1:9,100]
# REPEAT THE ABOVE TWO STEPS A FEW TIMES WITH DIFFERENT STARTING VALUES AND COMPARE THE PLOTS AND Lo
gLikelihood VALUES.
# AN EASY WAY TO GET DIFFERENT STARTING VALUES IS TO LET: aa_aa*(runif(9)+0.5)  
###################################################################################################


###################################################################################################
###################################################################################################