#############################################################
# Density function

dbbs = function(x, a=1, b=1, g=0){
  if(a<0){break} ;if(b<0){break}
  t=(sqrt(x/b)-sqrt(b/x))/a
  (x+b)*dnorm(abs(t)+g)/(4*a*sqrt(b)*(x^1.5)*pnorm(-g))
}

#############################################################
# Distribution function

pbbs = function(x,a=1,b=1,g=0){
  if(a<0){break} ;if(b<0){break}
  t=(sqrt(x/b)-sqrt(b/x))/a
  si = ifelse(x<b,1,0)
  (pnorm(t-g)/(2*pnorm(-g)))^si * 
    (.5 + (pnorm(t+g)-pnorm(g))/(2*pnorm(-g)))^(1-si)
}

#############################################################
# Random number generator
n = 10
i = 1

rbbs = function(n,a=1,b=1,g=0){
  if(a<0){break} ;if(b<0){break}
  u = runif(n,0,1)
  t = rep(0,n)
  for(i in 1:n){
    if(u[i]<0.5){t[i] = qnorm(2*pnorm(-g)*u[i]) + g}
    else{t[i] = qnorm(2*pnorm(-g)*(u[i]-1) + 1) - g} 
  }
  x = b*((t*a/2) + sqrt((t*a/2)^2 + 1))^2
  return(x)
}

#############################################################
# Moments of the BBS(a,b,g)

bbs_moments = function(r,a,b,g){
  if(a<0){break} ;if(b<0){break}
  E = 0
  for(k in 0:r){
    for(j in 0:k){
      m = 2*(r+j-k)
      for(s in 0:m){
        E = E + choose(2*r,2*k)*choose(k,j)*choose(m,s)*((a/2)^m)*
          ((-g)^(m-s))*dra(s,g) 
      }
    }
  }
  E = (b^r)*E/pnorm(-g)
  return(E)
}

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

# Incomplete moments of the standard normal distribution
d0 = function(a) 1 - pnorm(a)
d1 = function(a) dnorm(a)
d2 = function(a) 1 + a*dnorm(a) - pnorm(a)
d3 = function(a) (2+a*a)*dnorm(a)
d4 = function(a) 3+3*a*dnorm(a)+(a^3)*dnorm(a)-3*pnorm(a)
dra = function(r,a){
  f = function(x) (x^r)*dnorm(x)
  integrate(f,a,Inf)$value
}

# Second moment functions
k1 = function(a, g) ((a*g)^4)/2 + 2*(a*g)^2 + 1
k2 = function(a, g) 2*a*a*g*((a*g)^2 + 2)
k3 = function(a, g) a*a*(3*(a*g)^2 + 2)

# Function to estimate gamma in moment method
m2ag = function(a,g,u,r,v){
  (v - u*r*(k1(a,g)*d0(g) - k2(a,g)*d1(g) + k3(a,g)*d2(g) - 2*(a^4)*g*d3(g) + 
              (a^4)*d4(g)/2)/pnorm(-g))
}

#############################################################
# Moment method (MM) of estimation

# This function returns the estimatives of a BBS distribution fit for the
# data vector x
MMEbbs = function(x){
  n = length(x)
  u = mean(x)
  v = sum(x*x)/n
  r = n/sum(1/x)
  
  # Estimative of beta
  b = sqrt(u*r)
  
  g = 0 # initial guess of gamma
#  cont = 0; dif = 0
  for(i in 1:10){
    a = sqrt(2*(sqrt(u/r)-1)/(1+g*g-g*dnorm(g)/pnorm(g)))
#    if(sign(m2ag(a,-10,u,r,v))==sign(m2ag(a,10,u,r,v)) || cont>1){break}
    if(sign(m2ag(a,-10,u,r,v))==sign(m2ag(a,10,u,r,v))){break}
     res = uniroot(function(x) m2ag(a,x,u,r,v), c(-10,10),tol=.01,maxiter=1000)
     g = res$root
#    gn = res$root
#    if(abs(dif) < abs(gn-g)){cont=cont+1};dif = gn - g
#    g = gn
  }
  # The values of gamma and alpha are estimated iteratively. Gamma is obtained finding
  # the root of the function m2ag
  
  list(alpha = a, beta = b, gamma = g)
}

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

# t function
tf = function(x,a,b) (sqrt(x/b)-sqrt(b/x))/a
# w function
wf = function(g) dnorm(g)/pnorm(-g)
# Expectation E((X+beta)^(-2))
IK2 = function(a, b, g){
  fk2 = function(x){
    t=tf(x,a,b)
    dnorm(abs(t)+g)/(4*(x+b)*a*sqrt(b)*(x^1.5)*pnorm(-g))
  }
  integrate(fk2,0,Inf)$value
}

# Log-likelihood function
loglik = function(x, a, b, g) sum(log(dbbs(x,a,b,g)))

# tf(.5,1,1);wf(0);IK2(1,1,1) #tests

# Log-likelihood first derivatives
# alpha
dlda = function(x, a, b, g){
  t = tf(x,a,b)
  sum( -1 + (t^2) - g*abs(t) )/a
} 

# beta
dldb = function(x, a, b, g){
  t = tf(x,a,b)
  sum( -1/(2*b) + 1/(x+b) + sign(t)*(abs(t)+g)*(sqrt(x)+b/sqrt(x))/(2*a*(b^1.5)) )
} 

# gamma
dldg = function(x, a, b, g){
  t = tf(x,a,b)
  n = length(x)
  n*wf(g) - sum(abs(t)) - n*g
} 

# Expectation of the log-likelihood second derivatives
# n is the vector length


# alpha and alpha
d2lda2 = function(a, g, n){
  w = wf(g)
  -n*( 2 + g*(g-w) )/(a^2)
} 

# alpha and gamma
d2ldadg = function(a, g, n){
  w = wf(g)
  -n*(g-w)/a
} 

# beta and beta
d2ldb2 = function(a, b, g, n){
  w = wf(g)
  Kb = IK2(a, b, g)
  -n*( Kb + 1/(a*b)^2 + g*(g-w)/(4*(b^2)) )
}

# gamma and gamma
d2ldg2 = function(a, g, n){
  w = wf(g)
  -n*(w*(g-w) + 1)
}

# The remaining expectations are zero

# Function to calculate the Score vector
ScrVec = function(x,a,b,g) c(dlda(x,a,b,g), dldb(x,a,b,g), dldg(x,a,b,g))

# Observed Matrix

HessMat = function(x, a, b, g){
  t = (sqrt(x/b) - sqrt(b/x))/a
  w = wf(g); n = length(x)
  jaa = sum(1 - 3*t^2 - 2*g*abs(t))/a^2
  jab = sum(-((x/b^2)-(1/x))/a^3 - g*sign(t)*(sqrt(x)+(b/sqrt(x)))/(2*(a^2)*(b^1.5)))
  jbb = sum( 1/(2*b^2) - 1/(x+b)^2 - x/((a^2)*(b^3)) - 
               g*sign(t)*(3*sqrt(x)+b/sqrt(x))/(4*a*(b^2.5)) )
  jag = sum(abs(t))/a
  jbg = sum(sign(t)*(sqrt(x)+b/sqrt(x))/(2*a*(b^1.5)))
  jgg = -n*(w*(g-w) + 1)
  matrix(c(jaa,jab,jag,jab,jbb,jbg,jag,jbg,jgg),3,3)
}

# Function to calculate the Fisher Information Matrix
InfMat = function(a,b,g,n){
  (-1)*matrix(c(d2lda2(a, g, n), 0, d2ldadg(a, g, n), 0, d2ldb2(a, b, g, n), 0,
            d2ldadg(a, g, n), 0, d2ldg2(a, g, n)), 3, 3)
}

# Maximum Likelihood Estimation (MLE) method
require(MASS)
MLEbbs = function(x, n.iter=15){
  # initial guess
  theta = MMEbbs(x)
  theta = c(theta$alpha, theta$beta, theta$gamma)
  n = length(x)
  
#  K = InfMat(theta[1],theta[2],theta[3],n);K
#  Ki = solve(InfMat(theta[1],theta[2],theta[3],n));Ki
  J = -(1)*HessMat(x,theta[1],theta[2],theta[3]);J
  Ji = solve(J);Ji
  U = ScrVec(x,theta[1],theta[2],theta[3]);U
#  theta0 = theta + Ki%*%U;theta0
  theta0 = theta + Ji%*%U;theta0
  theta = theta0
  
  for(i in 1:n.iter) theta = theta + ginv(InfMat(theta[1],theta[2],theta[3],n))%*%
    ScrVec(x,theta[1],theta[2],theta[3])
  
  list(alpha = theta[1], beta = theta[2], gamma = theta[3], 
       U = ScrVec(x,theta[1],theta[2],theta[3]),
       K = InfMat(theta[1],theta[2],theta[3],n))
}

# Maximum likelihood using the optim function
MLEbbsOptim = function(dat,theta = NA){
  # initial guess
  if(is.na(theta[1])){
  theta = MMEbbs(dat)
  theta = c(theta$alpha, theta$beta, theta$gamma)}

  fn = function(x,theta){-loglik(x,theta[1],theta[2],theta[3])}
  #grad = function(x,theta){-ScrVec(x,theta[1],theta[2],theta[3])}
  
  res = optim(par=theta, fn, method="BFGS",hessian = TRUE, x=dat)

  theta = res$par

  list(alpha = theta[1], beta = theta[2], gamma = theta[3], 
       U = ScrVec(dat,theta[1],theta[2],theta[3]),
       K = InfMat(theta[1],theta[2],theta[3],length(dat)),
       Hessian = (-1)*res$hessian)
}

# Maximum likelihood using maxLik
require(maxLik)
MLEbbsMaxLik = function(x,theta = NA){
  # initial guess
  if(is.na(theta[1])){
    theta = MMEbbs(x)
    theta = c(theta$alpha, theta$beta, theta$gamma)}
  
  fn = function(theta){log(dbbs(x,theta[1],theta[2],theta[3]))}
  grad = function(theta){ScrVec(x,theta[1],theta[2],theta[3])}
#  hess = function(theta){InfMat(theta[1],theta[2],theta[3],length(x))}
#  hess = function(theta){HessMat(x,theta[1],theta[2],theta[3])}

  res = maxLik(fn, grad, hess=NULL, theta)
  a = res$estimate[1]; b = res$estimate[2]; g =  res$estimate[3]
  list(alpha = a, beta = b, gamma = g, 
       U = ScrVec(x,a,b,g),
       K = InfMat(a,b,g,length(x)),
       Hessian = res$hessian)
}

