Nothing
# derivatives and pdf functions for Standard Skewed Normal distribution
# will be merged with "code in the paper.txt"
####################### derivatives used in the model estimation ############
## goal: obtain the 1st and 2nd derivatives of the log-pdf of Y
## Y is SSN (with 0 mean, 1 var, and skewness \gamma)
## i.e. derivatives of log(g(y,gamma)) derivatives wrt y and gamma(skewness)
### derivatives of standardized SN with mean 0, variance 1 and shape parameter a
## input: argument y
## shape parameter a
## output: a list with 6 components (function values and 5 derivatives)
D.SN <- function(y,a){
xi<-sqrt(2/pi)*a/sqrt(1+a^2)
omega<-sqrt(1-xi^2)
x = omega*y+xi
# pre-requisite function values
# derivatives using SN when location = 0, scale = 1, shape = a
pdf.simple <- 2*dnorm(x)*pnorm(x*a)
D1.simple <- 1/pi*(-x*pi*pdf.simple + a*sqrt(2*pi)*dnorm(sqrt(1+a^2)*x))
D11.simple <- -(1 - x^2) * pdf.simple - a * x/pi * (2+a^2) * sqrt(2*pi) * dnorm(sqrt(1+a^2)*x)
D2.simple <- x/pi*sqrt(2*pi)*dnorm(sqrt(1+a^2)*x)
D22.simple <- -1*a*x^2*D2.simple
D12.simple <- 1/pi*sqrt(2*pi)*dnorm(sqrt(1+a^2)*x)*(1-(1+a^2)*x^2)
# derivatives of SSN with mean 0, variance 1 and shape a
sn <- omega*pdf.simple # function value
sn1 <- omega*D1.simple*omega # 1st derivative wrt. y
sn11 <- omega^3*D11.simple # 2nd derivative wrt. y
xi_a<- sqrt(2/pi)*(1+a^2)^(-1.5)
w_a<- -xi/omega*xi_a
x_a<- y*w_a + xi_a
#1st derivative wrt. a
sn2 <- w_a*pdf.simple+ omega*(D1.simple*x_a + D2.simple)
# mixture 2nd derivative wrt. y and a
sn12 <- 2*omega*w_a*D1.simple + omega^2*(D11.simple*x_a + D12.simple)
xi_aa = sqrt(2/pi)*(-3)*a*(1+a^2)^(-2.5)
w_aa = -(1-xi^2)^(-1.5)*xi_a*xi_a - xi/sqrt(1-xi^2)*xi_aa
x_aa = y*w_aa + xi_aa
# 2nd derivative wrt. a
sn22 <- w_aa*pdf.simple + 2*w_a*(D1.simple*x_a+D2.simple) +
omega*((D11.simple*x_a + D12.simple) *x_a +
D1.simple*x_aa + D12.simple*x_a + D22.simple)
return(list(sn = sn, sn1 = sn1, sn11 = sn11,
sn12 = sn12, sn2 = sn2, sn22 = sn22))}
#########################################################
### reparameterization between shape and skewness########
#########################################################
# transform shape to skewness
skewness.cp <- function(alpha) {
delta <- alpha/ sqrt( 1 + alpha^2)
temp <- delta * sqrt(2/pi)
(4 - pi)/2 * (temp)^3 / (1 - temp^2)^(3/2)}
# transfrom skewness to shape
shape.dp <- function(gamma) {
b <- sqrt(2/pi)
A <- sign(gamma) * (abs(2 * gamma/(4 - pi)))^(1/3)
delta <- A/(b * sqrt(1 + A^2))
lambda <- delta/sqrt(1 - delta^2)
return(lambda)}
#derivatives of gamma wrt shape(alpha)
D.gamma <- function(alpha){
gamma1 <- -3*pi*(pi-4)*alpha^2*sqrt(2)/ sqrt((pi + pi*alpha^2 - 2*alpha^2)^5)
gamma2 <- 3*(3*pi*alpha^2 - 6*alpha^2 - 2*pi) * pi *(pi - 4) *
alpha*sqrt(2) /sqrt((pi + pi*alpha^2 - 2*alpha^2)^7)
return(list(D1 = gamma1, D2 = gamma2))}
## pdf of sn(mean = 0, var = 1, skewness = gamma)
g <-function(y, gamma, log = FALSE) {
alpha <- shape.dp(gamma)
delta <- alpha/sqrt(1+alpha^2)
beta <- sqrt(1-2*delta^2/pi)
a <- beta*y+sqrt(1-beta^2)*sign(delta)
if (log == FALSE){
return(2*beta*dnorm(a)*pnorm(alpha*a))
}
if (log == TRUE){
return(log(2) + log(beta) + dnorm(a, log=T) + pnorm(alpha * a, log.p=T))
}
}
##log(g(y, gamma)) derivatives wrt y and gamma(skewness)
D.lg <- function(y, gamma){
fv <- g(y, gamma)
a <- shape.dp(gamma)
D.temp <- D.SN(y,a)
D1 <- D.temp$sn1
D11 <- D.temp$sn11
temp <- D.gamma(alpha = a)
temp1 <- temp[[1]]
temp2 <- temp[[2]]
D2 <- D.temp$sn2/temp1
D12 <- D.temp$sn12/temp1
D22 <- (D.temp$sn22 - temp2 * D2)/(temp1^2)
return(list(D1 = D1/fv, D2 = D2/fv,
D12 = D12/fv - 1/(fv^2)*D1*D2, D11 = D11/fv - 1/(fv^2)*D1^2,
D22 = D22/fv - 1/(fv^2)*D2^2))
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.