Nothing
optimfunct <- function (rhopar, v, verbose = FALSE){
vv <- v$bigG %*% c(rhopar[1], rhopar[1]^2) - v$litg
value <- sum(vv^2)
if (verbose) cat("function:", value, "lambda:", rhopar[1], "\n")
value
}
optimfunct_eff <- function (rhopar, v, vcmat, verbose = FALSE){
vv <- v$bigG %*% c(rhopar[1], rhopar[1]^2) - v$litg
value <- t(vv) %*% vcmat %*% vv
if (verbose) cat("function:", value, "lambda:", rhopar[1], "\n")
value
}
gg_het <- function(Ws, u, n) {
ub <- Ws %*% u
ubb <-Ws %*% ub
diag1 <- ub * u
diag2 <-ub^2
Wst <- t(Ws)
Diag1 <- Diagonal(n, as.vector(diag1))
tra <- sum(diag(Ws %*% Diag1 %*% Wst))
ubbub <- crossprod(ubb,ub)
first <- (ubbub - tra)
Diag2 <- Diagonal(n, as.vector(diag2))
tttt <- sum(diag(Ws %*% Diag2 %*% Wst))
ubbubb <- crossprod(ubb, ubb)
second <- (ubbubb - tttt)
uubb <- crossprod(u, ubb)
ubub <- crossprod(ub, ub)
third <- uubb + ubub
ububb <- crossprod(ub, ubb)
bigG <- matrix(0, 2, 2)
bigG[, 1] <- c(2*as.numeric(first), as.numeric(third))/n
bigG[, 2] <- -c(as.numeric(second), as.numeric(ububb))/n
diag3 <- u^2
Diag3 <- Diagonal(n, as.vector(diag3))
aa <- sum(diag(Ws%*%Diag3%*%Wst))
gamma1 <- (ubub - aa)
uub <- crossprod(u,ub)
litg <- c(as.numeric(gamma1), as.numeric(uub))/n
list(bigG = bigG, litg = litg)
}
gg_hom<-function (Ws, u, n) {
wu <- Ws %*% u
wwu<- Ws %*% wu
trwpw <- sum(Ws^2)
v.vec <- 1/ (1 + (trwpw/n)^2)
d <- (1/n)* trwpw
du <- d*u
ubub<- crossprod(wu)
udu<- crossprod(u,du)
uA1u<- v.vec * (ubub - udu)
uA2u<- crossprod(u, wu)
ubbubb <- crossprod(wwu)
dwu <- d*wu
ubdub <- crossprod(wu,dwu)
wuA1wu <- v.vec * (ubbubb - ubdub)
wuA2wu <- crossprod(wu, wwu)
udub <- crossprod(u,dwu)
ububb <- crossprod(wu,wwu)
upA1A1wu <- 2 * v.vec *(ububb - udub)
upA2A2wu <- crossprod(u,wwu) + ubub
bigG <- matrix(0, 2, 2)
bigG[, 1] <- c(as.numeric(upA1A1wu), as.numeric(upA2A2wu))/n
bigG[, 2] <- -c(as.numeric(wuA1wu), as.numeric(wuA2wu))/n
#print(bigG)
litg <- c(as.numeric(uA1u), as.numeric(uA2u))/n
# print(litg)
list(bigG = bigG, litg = litg, trwpw = trwpw, wu = wu, wwu = wwu, d = d, v.vec = v.vec)
}
psirhorho_hom <- function(rho, residuals, Hmat, Zmat, Ws, d, v.vec){
n <- length(residuals)
epsilon<- residuals - rho * Ws %*% residuals
mu3<- sum(epsilon^3) /n
mu4<- sum(epsilon^4) /n
sigma2 <- crossprod(epsilon)/n
Zstar <- Zmat - rho * Ws %*% Zmat
Qhz <- crossprod(Hmat,Zstar)/n
Qhh <- crossprod(Hmat)/n
Qhhi <- solve(Qhh)
sec.part <- t(Qhz) %*% Qhhi %*% Qhz
Pmat <- Qhhi %*% Qhz %*% solve(sec.part)
Tmat <- Hmat %*% Pmat
trA2A2 <- sum(Ws^2)
trA2A2I <- Diagonal(n) * (trA2A2/n)
A1 <- v.vec * (crossprod(Ws) - trA2A2I)
A2 <- Ws
A1pluA1 <- A1 + t(A1)
A2pluA2 <- A2 + t(A2)
A1pluA1eps <- A1pluA1 %*% epsilon
A2pluA2eps <- A2pluA2 %*% epsilon
ZA1pluA1eps <- crossprod(Zstar, A1pluA1eps)
ZA2pluA2eps <- crossprod(Zstar, A2pluA2eps)
alpha1 <- -1*ZA1pluA1eps/n
alpha2 <- -1*ZA2pluA2eps/n
a.vec1 <- Tmat %*% alpha1
a.vec2 <- Tmat %*% alpha2
tr22 <- sum(diag(A2pluA2 %*% A2pluA2))/(2*n)
# tr22<-(2*trA2A2 + 2* t(as.vector(t(Ws))) %*% as.vector(Ws)) /(2*n)
tr11 <- sum(diag(A1pluA1 %*% A1pluA1))/(2*n)
tr12 <- sum(diag(A1pluA1 %*% A2pluA2))/(2*n)
a1.a1 <- crossprod(a.vec1) /n
a2.a2 <- crossprod(a.vec2) /n
a1.a2 <- crossprod(a.vec1, a.vec2) /n
vecA1 <- diag(A1)
vecA2 <- diag(A2)
vecA1vecA1 <- crossprod(vecA1)/n
vecA2vecA2 <- crossprod(vecA2)/n
vecA1vecA2 <- crossprod(vecA1,vecA2)/n
a1.vecA1 <- crossprod(a.vec1, vecA1)
a1.vecA2 <- crossprod(a.vec1, vecA2)
a2.vecA2 <- crossprod(a.vec2, vecA2)
a2.vecA1 <- crossprod(a.vec2, vecA1)
la.term11 <- (a1.vecA1 + a1.vecA1)/n
la.term12 <- (a1.vecA2 + a2.vecA1)/n
la.term22 <- (a2.vecA2 + a2.vecA2)/n
effe <- mu4-3*sigma2^2
phi11 <- as.numeric(sigma2^2 * tr11 + sigma2 * a1.a1 + effe * vecA1vecA1 + mu3 * la.term11)
phi12 <- as.numeric(sigma2^2 * tr12 + sigma2 * a1.a2 + effe * vecA1vecA2 + mu3 * la.term12)
phi21 <- as.numeric(sigma2^2 * tr12 + sigma2 * a1.a2 + effe * vecA1vecA2 + mu3 * la.term12)
phi22 <- as.numeric(sigma2^2 * tr22 + sigma2 * a2.a2 + effe * vecA2vecA2 + mu3 * la.term22)
Phi <- cbind(rbind(phi11,phi21),rbind(phi12,phi22))
Phiinv <- solve(Phi)
list(Phi = Phi, Phiinv = Phiinv, Pmat= Pmat, A1 = A1, A2 =A2, a.vec1 = a.vec1, a.vec2 = a.vec2, epsilon = epsilon, Zstar = Zstar )
}
psirhorho_het<-function(rho, residuals, Hmat, Zmat, Ws, step1.c){
n <- length(residuals)
if(step1.c) Zstar <- Zmat
else Zstar<- Zmat - rho * Ws %*% Zmat
# print(dim(Zstar))
epsilon <- residuals - rho * Ws %*% residuals
Qhz <- crossprod(Hmat,Zstar)/n
Qhh <- crossprod(Hmat)/n
Qhhi <- solve(Qhh)
sec.part <- t(Qhz) %*% Qhhi %*% Qhz
Pmat <- Qhhi %*% Qhz %*% solve(sec.part)
Tmat <- Hmat %*% Pmat
trA2A2 <- sum(Ws^2)
A1 <- crossprod(Ws)
diag(A1) <- 0
A2 <- Ws
A1pluA1 <- A1 + t(A1)
A2pluA2 <- A2 + t(A2)
A1pluA1eps <- A1pluA1 %*% epsilon
A2pluA2eps <- A2pluA2 %*% epsilon
IrWp <- t(Diagonal(n) - rho * Ws)
ZpIrWp <- crossprod(Zmat,IrWp)
ZA1pluA1eps <- ZpIrWp %*% A1pluA1eps
ZA2pluA2eps <- ZpIrWp %*% A2pluA2eps
alpha1 <- -1*ZA1pluA1eps/n
alpha2 <- -1*ZA2pluA2eps/n
if(step1.c){
IrWpsi <- solve(Diagonal(n) - rho * t(Ws))
a.vec1 <- IrWpsi%*% Tmat %*% alpha1
a.vec2 <- IrWpsi%*% Tmat %*% alpha2
}
else{
a.vec1 <- Tmat %*% alpha1
a.vec2 <- Tmat %*% alpha2
}
gamma <- epsilon^2
gammas <- Diagonal(n, as.numeric(gamma))
#diag(gammas) <- gamma
a1Ga1 <- t(a.vec1) %*% gammas%*% a.vec1
a2Ga2 <- t(a.vec2) %*% gammas%*% a.vec2
a1Ga2 <- t(a.vec1) %*% gammas%*% a.vec2
tr1 <- sum(diag(A1pluA1%*%gammas%*%A1pluA1%*%gammas))/2
tr2 <- sum(diag(A2pluA2%*%gammas%*%A2pluA2%*%gammas))/2
tr3 <- sum(diag(A1pluA1%*%gammas%*%A2pluA2%*%gammas))/2
phi11 <- as.matrix((tr1 + a1Ga1)/n)
phi22 <- as.matrix((tr2 + a2Ga2)/n)
phi12 <- as.matrix((tr3 + a1Ga2)/n)
phir1 <- cbind(phi11, phi12)
phir2 <- cbind(phi12, phi22)
Phi <- rbind(as.matrix(phir1), as.matrix(phir2))
Phiinv <- solve(as.matrix(Phi))
list(Phi = Phi, Phiinv = Phiinv, Pmat= Pmat, A1 = A1, A2 =A2, a.vec1 = a.vec1, a.vec2 = a.vec2, epsilon = epsilon, Zstar = Zstar )
}
Omega_hom <- function(rho, Pmat, A1, A2, a.vec1, a.vec2, Hmat, bigG, Phirri, epsilon, Zstar){
n <- length(epsilon)
Jota <- bigG %*% c(1,2*rho)
mu3 <- sum(epsilon^3) / n
sigma2 <- crossprod(epsilon)/n
Qhh <- crossprod(Hmat)/n
a.vec <- cbind(as.numeric(a.vec1), as.numeric(a.vec2))
a.vec.d <- cbind(diag(A1), diag(A2))
Psidd <- as.numeric(sigma2) * Qhh
Psidr1<- as.numeric(sigma2)/n * crossprod(Hmat, a.vec)
Psidr2<- mu3/n * crossprod(Hmat, a.vec.d)
Psidr<- Psidr1 + Psidr2
Omegadd<- as.matrix(t(Pmat) %*% Psidd %*% Pmat)
Omegarr<- as.matrix(solve(t(Jota) %*% Phirri %*% Jota))
Omegadr<- as.matrix(t(Pmat) %*% Psidr %*% Phirri %*% Jota %*% Omegarr)
Omega<- cbind(rbind(Omegadd, t(Omegadr)), rbind(Omegadr, Omegarr))/n
list(Omega = Omega)
}
Omega_het <- function(rho, Pmat, A1, A2, a.vec1, a.vec2, Hmat, bigG, Phirri, epsilon, Zstar, Ws, step1.c){
n <- length(epsilon)
gamma <- epsilon^2
#gammas <- Matrix(diag(n), sparse = TRUE)
# diag(gammas) <- gamma
gammas <- Diagonal(n, as.numeric(gamma))
Jota<- bigG %*% c(1,2*rho)
a.vec <- cbind(as.numeric(a.vec1), as.numeric(a.vec2))
if(step1.c){
IrW <- Diagonal(n) - rho * Ws
IrWp <- t(IrW)
Psidd <- (t(Hmat) %*% IrW %*% gammas %*% IrWp %*% Hmat)/n
Psidr<- (t(Hmat) %*% IrW %*% gammas %*% a.vec)/n
}
else{
Psidd <- (t(Hmat) %*% gammas %*% Hmat)/n
Psidr<- (t(Hmat) %*% gammas %*% a.vec)/n
}
Omegadd<- as.matrix(t(Pmat) %*% Psidd %*% Pmat)
Omegarr<- as.matrix(solve(t(Jota) %*% Phirri %*% Jota))
Omegadr<- as.matrix(t(Pmat) %*% Psidr %*% Phirri %*% Jota %*% Omegarr)
Omega<- cbind(rbind(Omegadd, t(Omegadr)), rbind(Omegadr, Omegarr))/n
list(Omega = Omega)
}
psirhorho_het_mod<-function(rho, residuals, Hmat, Zmat, Ws, step1.c){
n<-length(residuals)
if(step1.c) Zstar <- Zmat
else Zstar<- Zmat - rho * Ws %*% Zmat
epsilon<- residuals - rho * Ws %*% residuals
Qhz<-crossprod(Hmat,Zstar)/n
Qhh<-crossprod(Hmat)/n
Qhhi<-solve(Qhh)
sec.part<- t(Qhz) %*% Qhhi %*% Qhz
Pmat<-Qhhi %*% Qhz %*% solve(sec.part)
Tmat<-Hmat %*% Pmat
trA2A2 <- sum(Ws^2)
A1 <- crossprod(Ws)
diag(A1) <- 0
A2 <- Ws
A1pluA1 <- A1 + t(A1)
A2pluA2 <- A2 + t(A2)
A1pluA1eps <- A1pluA1 %*% epsilon
A2pluA2eps <- A2pluA2 %*% epsilon
IrWp <- t(Diagonal(n) - rho * Ws)
ZpIrWp <- crossprod(Zmat,IrWp)
ZA1pluA1eps <- ZpIrWp %*% A1pluA1eps
ZA2pluA2eps <- ZpIrWp %*% A2pluA2eps
alpha1 <- -1*ZA1pluA1eps/n
alpha2 <- -1*ZA2pluA2eps/n
if(step1.c){
IrWpsi <- solve(Diagonal(n) - rho * t(Ws))
a.vec1 <- IrWpsi%*% Tmat %*% alpha1
a.vec2 <- IrWpsi%*% Tmat %*% alpha2
}
else{
a.vec1 <- Tmat %*% alpha1
a.vec2 <- Tmat %*% alpha2
}
gamma<-epsilon^2
# gammas<-Matrix(diag(n), sparse = TRUE)
# diag(gammas) <- gamma
gammas <- Diagonal(n, as.numeric(gamma))
a1Ga1<-t(a.vec1) %*% gammas%*% a.vec1
a2Ga2<-t(a.vec2) %*% gammas%*% a.vec2
a1Ga2<-t(a.vec1) %*% gammas%*% a.vec2
tr1<-sum(diag(A1pluA1%*%gammas%*%A1pluA1%*%gammas))/2
tr2<-sum(diag(A2pluA2%*%gammas%*%A2pluA2%*%gammas))/2
tr3<-sum(diag(A1pluA1%*%gammas%*%A2pluA2%*%gammas))/2
phi11<- as.matrix((tr1 )/n)
phi22<- as.matrix((tr2 )/n)
phi12<- as.matrix((tr3 )/n)
phir1<- cbind(phi11, phi12)
phir2<- cbind(phi12, phi22)
Phi<-rbind(as.matrix(phir1), as.matrix(phir2))
Phiinv<-solve(as.matrix(Phi))
list(Phi = Phi, Phiinv = Phiinv, Pmat= Pmat, A1 = A1, A2 =A2, a.vec1 = a.vec1, a.vec2 = a.vec2, epsilon = epsilon, Zstar = Zstar )
}
psirhorho_hom_mod <- function(rho, residuals, Hmat, Zmat, Ws, d, v.vec){
n <- length(residuals)
epsilon<- residuals - rho * Ws %*% residuals
mu3<- sum(epsilon^3) /n
mu4<- sum(epsilon^4) /n
sigma2<-crossprod(epsilon)/n
Zstar <- Zmat - rho * Ws %*% Zmat
Qhz<-crossprod(Hmat,Zstar)/n
Qhh<-crossprod(Hmat)/n
Qhhi<-solve(Qhh)
sec.part<- t(Qhz) %*% Qhhi %*% Qhz
Pmat<-Qhhi %*% Qhz %*% solve(sec.part)
Tmat<-Hmat %*% Pmat
trA2A2 <- sum(Ws^2)
trA2A2I<- Diagonal(n) * (trA2A2/n)
A1 <- v.vec * (crossprod(Ws) - trA2A2I)
A2 <- Ws
A1pluA1 <- A1 + t(A1)
A2pluA2 <- A2 + t(A2)
A1pluA1eps<- A1pluA1 %*% epsilon
A2pluA2eps<- A2pluA2 %*% epsilon
ZA1pluA1eps<- crossprod(Zstar, A1pluA1eps)
ZA2pluA2eps<- crossprod(Zstar, A2pluA2eps)
alpha1 <- -1*ZA1pluA1eps/n
alpha2 <- -1*ZA2pluA2eps/n
a.vec1 <- Tmat %*% alpha1
a.vec2 <- Tmat %*% alpha2
tr22<-sum(diag(A2pluA2 %*% A2pluA2))/(2*n)
# tr22<-(2*trA2A2 + 2* t(as.vector(t(Ws))) %*% as.vector(Ws)) /(2*n)
tr11<-sum(diag(A1pluA1 %*% A1pluA1))/(2*n)
tr12<-sum(diag(A1pluA1 %*% A2pluA2))/(2*n)
a1.a1 <-crossprod(a.vec1) /n
a2.a2 <-crossprod(a.vec2) /n
a1.a2 <-crossprod(a.vec1, a.vec2) /n
vecA1<- diag(A1)
vecA2<- diag(A2)
vecA1vecA1<- crossprod(vecA1)/n
vecA2vecA2<- crossprod(vecA2)/n
vecA1vecA2<- crossprod(vecA1,vecA2)/n
a1.vecA1<- crossprod(a.vec1, vecA1)
a1.vecA2<- crossprod(a.vec1, vecA2)
a2.vecA2<- crossprod(a.vec2, vecA2)
a2.vecA1<- crossprod(a.vec2, vecA1)
la.term11 <- (a1.vecA1 + a1.vecA1)/n
la.term12 <- (a1.vecA2 + a2.vecA1)/n
la.term22 <- (a2.vecA2 + a2.vecA2)/n
effe<- mu4-3*sigma2^2
# phi11<- as.numeric(sigma2^2 * tr11 )
# phi12<- as.numeric(sigma2^2 * tr12 )
# phi21<- as.numeric(sigma2^2 * tr12 )
# phi22<- as.numeric(sigma2^2 * tr22 )
phi11<- as.numeric(sigma2^2 * tr11 + effe * vecA1vecA1 )
phi12<- as.numeric(sigma2^2 * tr12 + effe * vecA1vecA2 )
phi21<- as.numeric(sigma2^2 * tr12 + effe * vecA1vecA2 )
phi22<- as.numeric(sigma2^2 * tr22 + effe * vecA2vecA2 )
Phi<-cbind(rbind(phi11,phi21),rbind(phi12,phi22))
Phiinv<-solve(Phi)
list(Phi = Phi, Phiinv = Phiinv, Pmat= Pmat, A1 = A1, A2 =A2, a.vec1 = a.vec1, a.vec2 = a.vec2, epsilon = epsilon, Zstar = Zstar )
}
Omega_hom_mod<-function(rho, Pmat, A1, A2, a.vec1, a.vec2, Hmat, bigG, Phirri, epsilon, Zstar){
n<-length(epsilon)
Jota<- bigG %*% c(1,2*rho)
mu3<-sum(epsilon^3) / n
sigma2<-crossprod(epsilon)/n
Qhh <- crossprod(Hmat)/n
a.vec <- cbind(as.numeric(a.vec1), as.numeric(a.vec2))
a.vec.d<-cbind(diag(A1), diag(A2))
Psidd <- as.numeric(sigma2) * Qhh
Psidr1<- as.numeric(sigma2)/n * crossprod(Hmat, a.vec)
Psidr2<- mu3/n * crossprod(Hmat, a.vec.d)
Psidr<- Psidr2
Omegadd<- as.matrix(t(Pmat) %*% Psidd %*% Pmat)
Omegarr<- as.matrix(solve(t(Jota) %*% Phirri %*% Jota))
Omegadr<- as.matrix(t(Pmat) %*% Psidr %*% Phirri %*% Jota %*% Omegarr)
Omega<- cbind(rbind(Omegadd, t(Omegadr)), rbind(Omegadr, Omegarr))/n
list(Omega = Omega)
}
Omega_het_mod<-function(rho, Pmat, A1, A2, a.vec1, a.vec2, Hmat, bigG, Phirri, epsilon, Zstar, Ws, step1.c){
n<-length(epsilon)
gamma<-epsilon^2
# gammas<-Matrix(diag(n), sparse = TRUE)
# diag(gammas) <- gamma
gammas <- Diagonal(n, as.numeric(gamma))
Jota<- bigG %*% c(1,2*rho)
a.vec <- cbind(as.numeric(a.vec1), as.numeric(a.vec2))
if(step1.c){
IrW <- Diagonal(n) - rho * Ws
IrWp <- t(IrW)
Psidd <- (t(Hmat) %*% IrW %*% gammas %*% IrWp %*% Hmat)/n
Psidr<- (t(Hmat) %*% IrW %*% gammas %*% a.vec)/n
}
else{
Psidd <- (t(Hmat) %*% gammas %*% Hmat)/n
Psidr<- (t(Hmat) %*% gammas %*% a.vec)/n
}
Omegadd<- as.matrix(t(Pmat) %*% Psidd %*% Pmat)
Omegarr<- as.matrix(solve(t(Jota) %*% Phirri %*% Jota))
Omegadr<- as.matrix(t(Pmat) %*% Psidr %*% Phirri %*% Jota %*% Omegarr)
# print(dim(Omegadr))
Omega<- cbind(rbind(Omegadd, matrix(0,1,2)), rbind(matrix(0,2,1), Omegarr))/n
list(Omega = Omega)
}
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.