The script you sent seems to be too large for me to spend time debugging,
plus it seems that there are files being read that are not part of the
data. I would suggest that if you have a long running program that you put
messages (cat('reached part 1\n')) throughout the program so you can trace
what it is doing. R typically does not stop on its own; it is only
following the directions of the program you have written and there is a bug
in it somewhere. So you will have to apply some basic debugging skills to
help you track it down.
There are a lot of things to look for. Does it take longer doing through
each loop (e.g., you have a vector that is growing and taking longer to
access). Is the system paging because you have run out of physical memory?
Do you have some condition in a 'while' loop that is never satisfied? Only
you will be able to determine this after you have put some trace statements
in your program and localized where the problem is.
On Sat, Mar 27, 2010 at 7:19 AM, Nomen Omen <h...@azet.sk> wrote:
> Dear addresses, I need perform a batch of 10 000 simulations for each of
> 4 options considered. (The idea is to obtain the parameter estimates in
> a heteroskedastic linear regression model - with additive or mixed
> heteroskedasticity - via the Kenward-Roger small-sample adjusted
> covariance matrix of disturbances). For this purpose I wrote an R
> program which would capture all possible options (true
> heteroskedasticity = additive or mixed / assumed heteroskedasticity =
> additive or mixed): simulate data, compute ML estimates of covariance
> components parameters (on the basis of maxLik package), and compute the
> resulting beta estimates. Simulations are performed by repeat command
> and, after every simulation step, the result is exported to an external
> csv file. Since sometimes an error was encountered during the simulation
> and the entire program broke down, I made use of the function TRY at
> critical evaluation steps to prevent the cessation of the program
> evaluation: R evaluates the result and if there should be a problem, it
> moves onto another step. Therefore, I distinguish between bad
> simulations and good simulations (and have two csv files for the
> results). There, however, another problem arises: R seems to be working
> (computing), it is not frozen or anything like that, and yet the saving
> to the files stops. I should say that it looks as if R was functioning
> fully but the simulations stopped. It is possible that R went into an
> undesired loop for no reasons that I know of. It requires of me to
> service all computers and, after the simulations have frozen, to restart
> the repeat command, which is not sound for I cannot be (or perhaps am
> not willing to be) with computers non-stop. The computers I use are of
> 2GB RAM and R is the version 2.10.1 (2009-12-14). I have no idea where
> the problem can lie and how to solve it. Can you please give me some
> suggestions how to handle this obstacle to speedy simulating? I shall be
> much thankful for any suggestion. Later I paste the example of the set
> of commands that I use, all though it might be not sensible, because it
> all is too complicated. I am convinced that of foremost importance is
> the last repeat command. Thank you very much. Faithfully yours, Martin
> Boïa (Matej Bel University / Faculty of Natural Sciences / Banská
> Bystrica / Slovakia).
>
> # zavádza sériu podporných funkcií
> # zavádza trace, sigma, fi, xs, xs, fixs, sxfixs
>
> tr <- function(M) sum(diag(M))
> sigma_ADD <- function(alpha) diag(as.vector(alpha %*% t(z)),
> nrow=nrow(z), ncol=nrow(z))
> sigma_MIX <- function(alpha) diag((as.vector(alpha %*% t(z)))^2,
> nrow=nrow(z), ncol=nrow(z))
> sigma_EXP <- function(alpha) diag(exp(as.vector(alpha %*% t(z))),
> nrow=nrow(z), ncol=nrow(z))
>
> fi <- function(sigma) solve(crossprod(x, crossprod(solve(sigma),x)))
>
> xs <- function(sigma) crossprod(x,solve(sigma))
> sx <- function(sigma) crossprod(solve(sigma),x)
> fixs <- function(sigma) crossprod(fi(sigma),xs(sigma))
> sxfixs <- function(sigma)
> crossprod(xs(sigma),crossprod(fi(sigma),xs(sigma)))
>
> # zavádza (ncol(z) x 1) vektor núl, ktorý má na i-tom mieste 1
>
> ei <- function(i) { ei = rep(0,ncol(z))
> ei[i] = 1
> return(ei) }
>
> # zavádza deriváciu dSIGMA/dalfa_i a zmie¹anú deriváciu
> (dSIGMA^2)/(ddalfa_i*dalfa_j)
> # zavádza maticový súèin (sigma)^(-1)*dSIGMA/dalfa_i
> # zavádza maticový súèin (sigma)^(-1)*dSIGMA/dalfa_i*(sigma)^(-1)
>
> dsigma.dalphai_ADD <- function(i) {
> diag(as.vector(ei(i) %*% t(z)), nrow=nrow(z), ncol=nrow(z)) }
> dsigma.dalphai_MIX <- function(sigma,i) {
> 2*crossprod(diag(as.vector(ei(i) %*% t(z)), nrow=nrow(z),
> ncol=nrow(z)),sigma^(1/2)) }
> dsigma.dalphai_EXP <- function(sigma,i) {
> crossprod(diag(as.vector(ei(i) %*% t(z)), nrow=nrow(z),
> ncol=nrow(z)),sigma) }
>
> d2sigma.dalphaij_MIX <- function(i,j) {
> eiz <- diag(as.vector(ei(i) %*% t(z)), nrow=nrow(z), ncol=nrow(z))
> ejz <- diag(as.vector(ei(j) %*% t(z)), nrow=nrow(z), ncol=nrow(z))
> d2sigma.dalphaij <- 2*crossprod(eiz,ejz)
> return(d2sigma.dalphaij) }
> d2sigma.dalphaij_EXP <- function(sigma,i,j) {
> eiz <- diag(as.vector(ei(i) %*% t(z)), nrow=nrow(z), ncol=nrow(z))
> ejz <- diag(as.vector(ei(j) %*% t(z)), nrow=nrow(z), ncol=nrow(z))
> d2sigma.dalphaij <- crossprod(crossprod(eiz,ejz),sigma)
> return(d2sigma.dalphaij) }
>
> sigmaDERIVi_ADD <- function(sigma,i)
> crossprod(solve(sigma),dsigma.dalphai_ADD(i))
> sigmaDERIVi_MIX <- function(sigma,i)
> crossprod(solve(sigma),dsigma.dalphai_MIX(sigma,i))
> sigmaDERIVi_EXP <- function(sigma,i)
> crossprod(solve(sigma),dsigma.dalphai_EXP(sigma,i))
>
> sigmaDERIVisigma_ADD <- function(sigma,i) {
> crossprod(solve(sigma),crossprod(dsigma.dalphai_ADD(i),solve(sigma))) }
> sigmaDERIVisigma_MIX <- function(sigma,i) {
> crossprod(solve(sigma),crossprod(dsigma.dalphai_MIX(sigma,i),solve(sigma)))
> }
> sigmaDERIVisigma_EXP <- function(sigma,i) {
> crossprod(solve(sigma),crossprod(dsigma.dalphai_EXP(sigma,i),solve(sigma)))
> }
>
>
> # zavádza pi, qij, rij (iba pre zmie¹ané a multiplikatívny model
> heteroskedasticity)
>
> pi_ADD <- function(sigma,i)
> -1*crossprod(x,crossprod(sigmaDERIVisigma_ADD(sigma,i),x))
> pi_MIX <- function(sigma,i)
> -1*crossprod(x,crossprod(sigmaDERIVisigma_MIX(sigma,i),x))
> pi_EXP <- function(sigma,i)
> -1*crossprod(x,crossprod(sigmaDERIVisigma_EXP(sigma,i),x))
>
> qij_ADD <- function(sigma,i,j) {
> SX <- sx(sigma)
> dSIGMA.dalphaj <- dsigma.dalphai_ADD(j)
> SIGMAderivi <- sigmaDERIVi_ADD(sigma,i)
> qij <- crossprod(SX,crossprod(SIGMAderivi,crossprod(dSIGMA.dalphaj,SX)))
> return(qij) }
> qij_MIX <- function(sigma,i,j) {
> SX <- sx(sigma)
> dSIGMA.dalphaj <- dsigma.dalphai_MIX(sigma,j)
> SIGMAderivi <- sigmaDERIVi_MIX(sigma,i)
> qij <- crossprod(SX,crossprod(SIGMAderivi,crossprod(dSIGMA.dalphaj,SX)))
> return(qij) }
> qij_EXP <- function(sigma,i,j) {
> SX <- sx(sigma)
> dSIGMA.dalphaj <- dsigma.dalphai_EXP(sigma,j)
> SIGMAderivi <- sigmaDERIVi_EXP(sigma,i)
> qij <- crossprod(SX,crossprod(SIGMAderivi,crossprod(dSIGMA.dalphaj,SX)))
> return(qij) }
>
> rij_MIX <- function(sigma,i,j)
> crossprod(sx(sigma),crossprod(d2sigma.dalphaij_MIX(i,j),sx(sigma)))
> rij_EXP <- function(sigma,i,j) {
> crossprod(sx(sigma),crossprod(d2sigma.dalphaij_EXP(sigma,i,j),sx(sigma)))
> }
>
> # zavádza funkciu pre stanovenie REML likelihood
> # je funkciou alpha conditional on x, y, z
>
>
> #########################################################################################
> reml.loglik_ADD <- function(alpha) { SIGMA <- sigma_ADD(alpha)
> a <- -0.5*log(det(SIGMA))
> b <- -0.5*log(det(crossprod(x,crossprod(solve(SIGMA),x))))
> c <-
> -0.5*as.numeric(crossprod(y,crossprod(solve(SIGMA)-sxfixs(SIGMA),y)))
> logl <- a + b + c
> return(logl) }
>
> reml.loglik_MIX <- function(alpha) { SIGMA <- sigma_MIX(alpha)
> a <- -0.5*log(det(SIGMA))
> b <- -0.5*log(det(crossprod(x,crossprod(solve(SIGMA),x))))
> c <-
> -0.5*as.numeric(crossprod(y,crossprod(solve(SIGMA)-sxfixs(SIGMA),y)))
> logl <- a + b + c
> return(logl) }
>
> reml.loglik_EXP <- function(alpha) { SIGMA <- sigma_EXP(alpha)
> a <- -0.5*log(det(SIGMA))
> b <- -0.5*log(det(crossprod(x,crossprod(solve(SIGMA),x))))
> c <-
> -0.5*as.numeric(crossprod(y,crossprod(solve(SIGMA)-sxfixs(SIGMA),y)))
> logl <- a + b + c
> return(logl) }
>
> #########################################################################################
>
> # zavádza funkciu pre stanovenie gradientného vektora pre REML
> likelihood
> # je funkciou alpha conditional on x, y, z
> # obsahuje podfunkciu dl.dalphai, ktorá spoèítava deriváciu
> likelihood function podµa alpha i
>
>
> #########################################################################################
> ######## gradientný vektor aditívneho modelu
> ############################################
>
> gradvec_ADD = function(alpha) { dl.dalphai <- function(i) {
>
> SIGMA <- sigma_ADD(alpha)
> FI <- fi(SIGMA)
> FIxs <- fixs(SIGMA)
> sxFIxs <- sxfixs(SIGMA)
> Pi <- pi_ADD(SIGMA,i)
>
> dSIGMA.dalphai <- dsigma.dalphai_ADD(i)
> SIGMAderiviSIGMA <- sigmaDERIVisigma_ADD(SIGMA,i)
>
> dl.dalphai.1 <- -1/2*tr(crossprod(solve(SIGMA),dSIGMA.dalphai ))
> dl.dalphai.2 <- -1/2*tr(crossprod(FI,Pi))
> dl.dalphai.3 <-
> +1/2*as.numeric(crossprod(y,crossprod(SIGMAderiviSIGMA,y)))
> INSIDE1 <- crossprod(FIxs,crossprod(Pi,FIxs))
> dl.dalphai.4 <- -1/2*as.numeric(crossprod(y,crossprod(INSIDE1,y)))
> INSIDE2 <- crossprod(solve(SIGMA),crossprod(dSIGMA.dalphai ,sxFIxs))
> dl.dalphai.5 <-
> +1/2*(-2)*as.numeric(crossprod(y,crossprod(t(INSIDE2),y)))
>
> dl.dalphai = dl.dalphai.1 + dl.dalphai.2 + dl.dalphai.3 + dl.dalphai.4 +
> dl.dalphai.5
>
> return(dl.dalphai)}
>
>
> grad = rep(dl.dalphai(1),c(1))
>
> for(i in 2:length(alpha)) {
> grad = append(grad,rep(dl.dalphai(i),c(1)), after = c(i)) }
>
> return(grad) }
>
> ######## gradientný vektor zmie¹aného modelu
> ############################################
>
> gradvec_MIX = function(alpha) { dl.dalphai <- function(i) {
>
> SIGMA <- sigma_MIX(alpha)
> FI <- fi(SIGMA)
> FIxs <- fixs(SIGMA)
> sxFIxs <- sxfixs(SIGMA)
> Pi <- pi_MIX(SIGMA,i)
>
> dSIGMA.dalphai <- dsigma.dalphai_MIX(SIGMA,i)
> SIGMAderiviSIGMA <- sigmaDERIVisigma_MIX(SIGMA,i)
>
> dl.dalphai.1 <- -1/2*tr(crossprod(solve(SIGMA),dSIGMA.dalphai ))
> dl.dalphai.2 <- -1/2*tr(crossprod(FI,Pi))
> dl.dalphai.3 <-
> +1/2*as.numeric(crossprod(y,crossprod(SIGMAderiviSIGMA,y)))
> INSIDE1 <- crossprod(FIxs,crossprod(Pi,FIxs))
> dl.dalphai.4 <- -1/2*as.numeric(crossprod(y,crossprod(INSIDE1,y)))
> INSIDE2 <- crossprod(solve(SIGMA),crossprod(dSIGMA.dalphai ,sxFIxs))
> dl.dalphai.5 <-
> +1/2*(-2)*as.numeric(crossprod(y,crossprod(t(INSIDE2),y)))
>
> dl.dalphai = dl.dalphai.1 + dl.dalphai.2 + dl.dalphai.3 + dl.dalphai.4 +
> dl.dalphai.5
>
> return(dl.dalphai)}
>
>
> grad = rep(dl.dalphai(1),c(1))
>
> for(i in 2:length(alpha)) {
> grad = append(grad,rep(dl.dalphai(i),c(1)), after = c(i)) }
>
> return(grad) }
>
> ######## gradientný vektor multiplikatívneho modelu
> #####################################
>
> gradvec_EXP = function(alpha) { dl.dalphai <- function(i) {
>
> SIGMA <- sigma_EXP(alpha)
> FI <- fi(SIGMA)
> FIxs <- fixs(SIGMA)
> sxFIxs <- sxfixs(SIGMA)
> Pi <- pi_EXP(SIGMA,i)
>
> dSIGMA.dalphai <- dsigma.dalphai_EXP(SIGMA,i)
> SIGMAderiviSIGMA <- sigmaDERIVisigma_EXP(SIGMA,i)
>
> dl.dalphai.1 <- -1/2*tr(crossprod(solve(SIGMA),dSIGMA.dalphai ))
> dl.dalphai.2 <- -1/2*tr(crossprod(FI,Pi))
> dl.dalphai.3 <-
> +1/2*as.numeric(crossprod(y,crossprod(SIGMAderiviSIGMA,y)))
> INSIDE1 <- crossprod(FIxs,crossprod(Pi,FIxs))
> dl.dalphai.4 <- -1/2*as.numeric(crossprod(y,crossprod(INSIDE1,y)))
> INSIDE2 <- crossprod(solve(SIGMA),crossprod(dSIGMA.dalphai ,sxFIxs))
> dl.dalphai.5 <-
> +1/2*(-2)*as.numeric(crossprod(y,crossprod(t(INSIDE2),y)))
>
> dl.dalphai = dl.dalphai.1 + dl.dalphai.2 + dl.dalphai.3 + dl.dalphai.4 +
> dl.dalphai.5
>
> return(dl.dalphai)}
>
>
> grad = rep(dl.dalphai(1),c(1))
>
> for(i in 2:length(alpha)) {
> grad = append(grad,rep(dl.dalphai(i),c(1)), after = c(i)) }
>
> return(grad) }
>
>
> #########################################################################################
>
> # zavádza funkciu pre stanovenie Hessovej matice pre REML likelihood
> # je funkciou alpha conditional on x, y, z
> # obsahuje podfunkciu d2l.dalphaij, èo spoèítava druhú deriváciu
> loglikelihood podµa alpha i a alpha j
>
>
> #########################################################################################
> ######## Hessova matica aditívneho modelu
> ##############################################
>
> hessmat_ADD = function(alpha) {
>
> hess = matrix(nrow = ncol(z), ncol = ncol(z))
>
> d2l.dalphaij <- function(i,j) {
>
> SIGMA <- sigma_ADD(alpha)
> FI <- fi(SIGMA)
> XS <- xs(SIGMA)
> SX <- sx(SIGMA)
> FIxs <- fixs(SIGMA)
> sxFI <- t(FIxs)
> sxFIxs <- sxfixs(SIGMA)
> Pi <- pi_ADD(SIGMA,i)
> Pj <- pi_ADD(SIGMA,j)
> Qij <- qij_ADD(SIGMA,i,j)
> dSIGMA.dalphai <- dsigma.dalphai_ADD(i)
> dSIGMA.dalphaj <- dsigma.dalphai_ADD(j)
> SIGMAderiviSIGMA <- sigmaDERIVisigma_ADD(SIGMA,i)
> SIGMAderivjSIGMA <- sigmaDERIVisigma_ADD(SIGMA,j)
> SIGMAderivi <- sigmaDERIVi_ADD(SIGMA,i)
> SIGMAderivj <- sigmaDERIVi_ADD(SIGMA,j)
>
> d2l.dalphaij.11 = 1/2*tr(crossprod(t(SIGMAderivi),SIGMAderivj))
> d2l.dalphaij.12 = 0
> d2l.dalpahij.21 = 1/2*tr(crossprod(FI,crossprod(Pi,crossprod(FI,Pj))))
> d2l.dalphaij.22 = -1*tr(crossprod(FI,Qij))
> d2l.dalphaij.23 = 0
> d2l.dalphaij.31 =
> -1*crossprod(y,crossprod(t(SIGMAderivj),crossprod(SIGMAderiviSIGMA,y)))
> d2l.dalphaij.32 = 0
> d2l.dalphaij.41 =
>
> crossprod(y,crossprod(t(SIGMAderivj),crossprod(FIxs,crossprod(Pi,crossprod(sxFI,y)))))
> d2l.dalphaij.42 =
>
> crossprod(y,crossprod(FIxs,crossprod(Pi,crossprod(FI,crossprod(Pj,crossprod(sxFI,y))))))
> d2l.dalphaij.43 =
> -1*crossprod(y,crossprod(FIxs,crossprod(Qij,crossprod(sxFI,y))))
> d2l.dalphaij.44 = 0
> d2l.dalphaij.51 =
>
> crossprod(y,crossprod(t(SIGMAderivj),crossprod(t(SIGMAderivi),crossprod(sxFIxs,y))))
> d2l.dalphaij.52 = 0
> d2l.dalphaij.53 =
>
> crossprod(y,crossprod(t(SIGMAderivi),crossprod(t(SIGMAderivj),crossprod(sxFIxs,y))))
> d2l.dalphaij.54 =
>
> crossprod(y,crossprod(t(SIGMAderivi),crossprod(FIxs,crossprod(Pj,crossprod(sxFI,y)))))
> d2l.dalphaij.55 =
>
> crossprod(y,crossprod(t(SIGMAderivi),crossprod(sxFIxs,crossprod(SIGMAderivj,y))))
>
> d2l.dalphaij = d2l.dalphaij.11 + d2l.dalphaij.12 + d2l.dalpahij.21 +
> d2l.dalphaij.22 + d2l.dalphaij.23 +
> d2l.dalphaij.31 + d2l.dalphaij.32 + d2l.dalphaij.41 + d2l.dalphaij.42 +
> d2l.dalphaij.43 +
> d2l.dalphaij.44 + d2l.dalphaij.51 + d2l.dalphaij.52 + d2l.dalphaij.53 +
> d2l.dalphaij.54 +
> d2l.dalphaij.55
>
> return(d2l.dalphaij)}
>
> for(i in 1:ncol(z)) { for(j in 1:ncol(z)) hess[i,j] = d2l.dalphaij(i,j)
> }
>
> return(hess)}
>
> ######## Hessova matica zmie¹aného modelu
> ##############################################
>
> hessmat_MIX = function(alpha) {
>
> hess = matrix(nrow = ncol(z), ncol = ncol(z))
>
> d2l.dalphaij <- function(i,j) {
>
> SIGMA <- sigma_MIX(alpha)
> FI <- fi(SIGMA)
> XS <- xs(SIGMA)
> SX <- sx(SIGMA)
> FIxs <- fixs(SIGMA)
> sxFI <- t(FIxs)
> sxFIxs <- sxfixs(SIGMA)
> Pi <- pi_MIX(SIGMA,i)
> Pj <- pi_MIX(SIGMA,j)
> Qij <- qij_MIX(SIGMA,i,j)
> Rij<- rij_MIX(SIGMA,i,j)
> dSIGMA.dalphai <- dsigma.dalphai_MIX(SIGMA,i)
> dSIGMA.dalphaj <- dsigma.dalphai_MIX(SIGMA,j)
> d2SIGMA.dalphaij <- d2sigma.dalphaij_MIX(i,j)
> SIGMAderiviSIGMA <- sigmaDERIVisigma_MIX(SIGMA,i)
> SIGMAderivjSIGMA <- sigmaDERIVisigma_MIX(SIGMA,j)
> SIGMAderivi <- sigmaDERIVi_MIX(SIGMA,i)
> SIGMAderivj <- sigmaDERIVi_MIX(SIGMA,j)
>
> d2l.dalphaij.11 = 1/2*tr(crossprod(t(SIGMAderivi),SIGMAderivj))
> d2l.dalphaij.12 = -1/2*tr(crossprod(solve(SIGMA),d2SIGMA.dalphaij))
> d2l.dalpahij.21 = 1/2*tr(crossprod(FI,crossprod(Pi,crossprod(FI,Pj))))
> d2l.dalphaij.22 = -1*tr(crossprod(FI,Qij))
> d2l.dalphaij.23 = 1/2*tr(crossprod(FI,Rij))
> d2l.dalphaij.31 =
> -1*crossprod(y,crossprod(t(SIGMAderivj),crossprod(SIGMAderiviSIGMA,y)))
> d2l.dalphaij.32 =
>
> 1/2*crossprod(y,crossprod(solve(SIGMA),crossprod(d2SIGMA.dalphaij,crossprod(solve(SIGMA),y))))
> d2l.dalphaij.41 =
>
> crossprod(y,crossprod(t(SIGMAderivj),crossprod(FIxs,crossprod(Pi,crossprod(sxFI,y)))))
> d2l.dalphaij.42 =
>
> crossprod(y,crossprod(FIxs,crossprod(Pi,crossprod(FI,crossprod(Pj,crossprod(sxFI,y))))))
> d2l.dalphaij.43 =
> -1*crossprod(y,crossprod(FIxs,crossprod(Qij,crossprod(sxFI,y))))
> d2l.dalphaij.44 =
> 1/2*crossprod(y,crossprod(FIxs,crossprod(Rij,crossprod(sxFI,y))))
> d2l.dalphaij.51 =
>
> crossprod(y,crossprod(t(SIGMAderivj),crossprod(t(SIGMAderivi),crossprod(sxFIxs,y))))
> d2l.dalphaij.52 =
>
> -1*crossprod(y,crossprod(solve(SIGMA),crossprod(d2SIGMA.dalphaij,crossprod(sxFIxs,y))))
> d2l.dalphaij.53 =
>
> crossprod(y,crossprod(t(SIGMAderivi),crossprod(t(SIGMAderivj),crossprod(sxFIxs,y))))
> d2l.dalphaij.54 =
>
> crossprod(y,crossprod(t(SIGMAderivi),crossprod(FIxs,crossprod(Pj,crossprod(sxFI,y)))))
> d2l.dalphaij.55 =
>
> crossprod(y,crossprod(t(SIGMAderivi),crossprod(sxFIxs,crossprod(SIGMAderivj,y))))
>
> d2l.dalphaij = d2l.dalphaij.11 + d2l.dalphaij.12 + d2l.dalpahij.21 +
> d2l.dalphaij.22 + d2l.dalphaij.23 +
> d2l.dalphaij.31 + d2l.dalphaij.32 + d2l.dalphaij.41 + d2l.dalphaij.42 +
> d2l.dalphaij.43 +
> d2l.dalphaij.44 + d2l.dalphaij.51 + d2l.dalphaij.52 + d2l.dalphaij.53 +
> d2l.dalphaij.54 +
> d2l.dalphaij.55
>
> return(d2l.dalphaij)}
>
> for(i in 1:ncol(z)) { for(j in 1:ncol(z)) hess[i,j] = d2l.dalphaij(i,j)
> }
>
> return(hess)}
>
> ######## Hessova matica multiplikatívneho modelu
> #######################################
>
> hessmat_EXP = function(alpha) {
>
> hess = matrix(nrow = ncol(z), ncol = ncol(z))
>
> d2l.dalphaij <- function(i,j) {
>
> SIGMA <- sigma_EXP(alpha)
> FI <- fi(SIGMA)
> XS <- xs(SIGMA)
> SX <- sx(SIGMA)
> FIxs <- fixs(SIGMA)
> sxFI <- t(FIxs)
> sxFIxs <- sxfixs(SIGMA)
> Pi <- pi_EXP(SIGMA,i)
> Pj <- pi_EXP(SIGMA,j)
> Qij <- qij_EXP(SIGMA,i,j)
> Rij<- rij_EXP(SIGMA,i,j)
> dSIGMA.dalphai <- dsigma.dalphai_EXP(SIGMA,i)
> dSIGMA.dalphaj <- dsigma.dalphai_EXP(SIGMA,j)
> d2SIGMA.dalphaij <- d2sigma.dalphaij_EXP(SIGMA,i,j)
> SIGMAderiviSIGMA <- sigmaDERIVisigma_EXP(SIGMA,i)
> SIGMAderivjSIGMA <- sigmaDERIVisigma_EXP(SIGMA,j)
> SIGMAderivi <- sigmaDERIVi_EXP(SIGMA,i)
> SIGMAderivj <- sigmaDERIVi_EXP(SIGMA,j)
>
>
>
> d2l.dalphaij.11 = 1/2*tr(crossprod(t(SIGMAderivi),SIGMAderivj))
> d2l.dalphaij.12 = -1/2*tr(crossprod(solve(SIGMA),d2SIGMA.dalphaij))
> d2l.dalpahij.21 = 1/2*tr(crossprod(FI,crossprod(Pi,crossprod(FI,Pj))))
> d2l.dalphaij.22 = -1*tr(crossprod(FI,Qij))
> d2l.dalphaij.23 = 1/2*tr(crossprod(FI,Rij))
> d2l.dalphaij.31 =
> -1*crossprod(y,crossprod(t(SIGMAderivj),crossprod(SIGMAderiviSIGMA,y)))
> d2l.dalphaij.32 =
>
> 1/2*crossprod(y,crossprod(solve(SIGMA),crossprod(d2SIGMA.dalphaij,crossprod(solve(SIGMA),y))))
> d2l.dalphaij.41 =
>
> crossprod(y,crossprod(t(SIGMAderivj),crossprod(FIxs,crossprod(Pi,crossprod(sxFI,y)))))
> d2l.dalphaij.42 =
>
> crossprod(y,crossprod(FIxs,crossprod(Pi,crossprod(FI,crossprod(Pj,crossprod(sxFI,y))))))
> d2l.dalphaij.43 =
> -1*crossprod(y,crossprod(FIxs,crossprod(Qij,crossprod(sxFI,y))))
> d2l.dalphaij.44 =
> 1/2*crossprod(y,crossprod(FIxs,crossprod(Rij,crossprod(sxFI,y))))
> d2l.dalphaij.51 =
>
> crossprod(y,crossprod(t(SIGMAderivj),crossprod(t(SIGMAderivi),crossprod(sxFIxs,y))))
> d2l.dalphaij.52 =
>
> -1*crossprod(y,crossprod(solve(SIGMA),crossprod(d2SIGMA.dalphaij,crossprod(sxFIxs,y))))
> d2l.dalphaij.53 =
>
> crossprod(y,crossprod(t(SIGMAderivi),crossprod(t(SIGMAderivj),crossprod(sxFIxs,y))))
> d2l.dalphaij.54 =
>
> crossprod(y,crossprod(t(SIGMAderivi),crossprod(FIxs,crossprod(Pj,crossprod(sxFI,y)))))
> d2l.dalphaij.55 =
>
> crossprod(y,crossprod(t(SIGMAderivi),crossprod(sxFIxs,crossprod(SIGMAderivj,y))))
>
> d2l.dalphaij = d2l.dalphaij.11 + d2l.dalphaij.12 + d2l.dalpahij.21 +
> d2l.dalphaij.22 + d2l.dalphaij.23 +
> d2l.dalphaij.31 + d2l.dalphaij.32 + d2l.dalphaij.41 + d2l.dalphaij.42 +
> d2l.dalphaij.43 +
> d2l.dalphaij.44 + d2l.dalphaij.51 + d2l.dalphaij.52 + d2l.dalphaij.53 +
> d2l.dalphaij.54 +
> d2l.dalphaij.55
>
> return(d2l.dalphaij)}
>
> for(i in 1:ncol(z)) { for(j in 1:ncol(z)) hess[i,j] = d2l.dalphaij(i,j)
> }
>
> return(hess)}
>
>
> #########################################################################################
> ######## zavádza information matrices a kovarianèné matice
> aditívneho modelu ############
>
> # expected information matrix
>
> expmat_ADD <- function(alpha) { mat = matrix(nrow = ncol(z), ncol =
> ncol(z))
>
> matij <- function(i,j) {
>
> SIGMA <- sigma_ADD(alpha)
> FI <- fi(SIGMA)
> Pi <- pi_ADD(SIGMA,i)
> Pj <- pi_ADD(SIGMA,j)
> Qij <- qij_ADD(SIGMA,i,j)
> SIGMAderivi <- sigmaDERIVi_ADD(SIGMA,i)
> SIGMAderivj <- sigmaDERIVi_ADD(SIGMA,j)
>
> matij.1 = 1/2*tr(crossprod(t(SIGMAderivi),SIGMAderivj))
> matij.2 =
> -1*tr(crossprod(FI,Qij)-1/2*crossprod(FI,crossprod(Pi,crossprod(FI,Pj))))
>
> matij = matij.1 + matij.2
>
> return(matij)}
>
> for(i in 1:ncol(z)) { for(j in 1:ncol(z)) mat[i,j] = matij(i,j) }
>
> return(mat)}
>
> # observed information matrix & average infromation matrix
>
> obsmat_ADD <- function(alpha) -1*hessmat_ADD(alpha)
> avgmat_ADD <- function(alpha) 1/2*(obsmat_ADD(alpha) +
> expmat_ADD(alpha))
>
> expcovmat_ADD <- function(alpha) solve(expmat_ADD(alpha))
> obscovmat_ADD <- function(alpha) solve(obsmat_ADD(alpha))
> avgcovmat_ADD <- function(alpha) solve(avgmat_ADD(alpha))
>
> ######## zavádza information matrices a kovarianèné matice
> zmie¹aného modelu ############
>
> # expected information matrix
>
> expmat_MIX = function(alpha) { mat = matrix(nrow = ncol(z), ncol =
> ncol(z))
>
> matij <- function(i,j) {
>
> SIGMA <- sigma_MIX(alpha)
> FI <- fi(SIGMA)
> Pi <- pi_MIX(SIGMA,i)
> Pj <- pi_MIX(SIGMA,j)
> Qij <- qij_MIX(SIGMA,i,j)
> SIGMAderivi <- sigmaDERIVi_MIX(SIGMA,i)
> SIGMAderivj <- sigmaDERIVi_MIX(SIGMA,j)
>
> matij.1 = 1/2*tr(crossprod(t(SIGMAderivi),SIGMAderivj))
> matij.2 =
> -1*tr(crossprod(FI,Qij)-1/2*crossprod(FI,crossprod(Pi,crossprod(FI,Pj))))
>
> matij = matij.1 + matij.2
>
> return(matij)}
>
> for(i in 1:ncol(z)) { for(j in 1:ncol(z)) mat[i,j] = matij(i,j) }
>
> return(mat)}
>
> # observed information matrix & average infromation matrix
>
> obsmat_MIX <- function(alpha) -1*hessmat_MIX(alpha)
> avgmat_MIX <- function(alpha) 1/2*(obsmat_MIX(alpha) +
> expmat_MIX(alpha))
>
> expcovmat_MIX <- function(alpha) solve(expmat_MIX(alpha))
> obscovmat_MIX <- function(alpha) solve(obsmat_MIX(alpha))
> avgcovmat_MIX <- function(alpha) solve(avgmat_MIX(alpha))
>
> ######## zavádza information matrices a kovarianèné matice
> mutliplikatívneho modelu #####
>
> # expected information matrix
>
> expmat_EXP = function(alpha) { mat = matrix(nrow = ncol(z), ncol =
> ncol(z))
>
> matij <- function(i,j) {
>
> SIGMA <- sigma_EXP(alpha)
> FI <- fi(SIGMA)
> Pi <- pi_EXP(SIGMA,i)
> Pj <- pi_EXP(SIGMA,j)
> Qij <- qij_EXP(SIGMA,i,j)
> SIGMAderivi <- sigmaDERIVi_EXP(SIGMA,i)
> SIGMAderivj <- sigmaDERIVi_EXP(SIGMA,j)
>
> matij.1 = 1/2*tr(crossprod(t(SIGMAderivi),SIGMAderivj))
> matij.2 =
> -1*tr(crossprod(FI,Qij)-1/2*crossprod(FI,crossprod(Pi,crossprod(FI,Pj))))
>
> matij = matij.1 + matij.2
>
> return(matij)}
>
> for(i in 1:ncol(z)) { for(j in 1:ncol(z)) mat[i,j] = matij(i,j) }
>
> return(mat)}
>
> # observed information matrix & average infromation matrix
>
> obsmat_EXP <- function(alpha) -1*hessmat_EXP(alpha)
> avgmat_EXP <- function(alpha) 1/2*(obsmat_EXP(alpha) +
> expmat_EXP(alpha))
>
> expcovmat_EXP <- function(alpha) solve(expmat_EXP(alpha))
> obscovmat_EXP <- function(alpha) solve(obsmat_EXP(alpha))
> avgcovmat_EXP <- function(alpha) solve(avgmat_EXP(alpha))
>
>
> #########################################################################################
>
> # zavádza kni¾nicu maxLik
>
> library(maxLik)
>
> # zavedie ¹tartovacie hodnoty jednotky pre alpha
>
> beta_ols <- function(x,y) crossprod(solve(crossprod(x)), crossprod(x,y))
> resids_ols <- function(x,y) (y - crossprod(t(x),beta_ols(x,y)))
> alpha_ols_ADD <- function(x,y)
> crossprod(solve(crossprod(z)),crossprod(z,resids_ols(x,y)^2))
> alpha_ols_MIX <- function(x,y)
> crossprod(solve(crossprod(z)),crossprod(z,abs(resids_ols(x,y))))
> alpha_ols_EXP <- function(x,y)
> crossprod(solve(crossprod(z)),crossprod(z,log(resids_ols(x,y)^2)))
> coefolsEST_ADD <- function(x,y) t(alpha_ols_ADD(x,y))
> coefolsEST_MIX <- function(x,y) t(alpha_ols_MIX(x,y))
> coefolsEST_EXP <- function(x,y) t(alpha_ols_EXP(x,y))
>
>
> kenward_roger_both <- function(model,estimated_alpha,beta_tested,x,y,z)
> { x = x; y = y; z = z
> if (model == "ADD") {
> inner <- function(i,j,alpha) {
> FI <- fi(sigma_ADD(alpha))
> Pi <- pi_ADD(sigma_ADD(alpha),i)
> Pj <- pi_ADD(sigma_ADD(alpha),j)
> Qij <- qij_ADD(sigma_ADD(alpha),i,j)
> W <- expcovmat_ADD(alpha)
> inner = matrix(nrow = ncol(x), ncol = ncol(x))
> inner = W[i,j]*(Qij-crossprod(Pi,crossprod(FI,Pj)))
> return(inner) }
> INNER = matrix(0, nrow = ncol(x), ncol = ncol(x))
> for(i in 1:ncol(z)) { for(j in 1:ncol(z)) INNER = INNER +
> inner(i,j,estimated_alpha) }
> a1a2ij <- function(i,j,alpha) {
> FI <- fi(sigma_ADD(alpha))
> Pi <- pi_ADD(sigma_ADD(alpha),i)
> Pj <- pi_ADD(sigma_ADD(alpha),j)
> W <- expcovmat_ADD(alpha)
> a1a2 = vector(mode = "numeric",length = 2)
> a1a2[1] <- W[i,j]*tr(crossprod(FI,Pi))*tr(crossprod(FI,Pj))
> a1a2[2] <- W[i,j]*tr(crossprod(crossprod(Pi,FI),crossprod(FI,Pj)))
> return(a1a2) }
> A1A2 = vector(mode = "numeric", length = 2)
> for(i in 1:ncol(z)) { for(j in 1:ncol(z)) A1A2 = A1A2 +
> a1a2ij(i,j,estimated_alpha) }
> A1 = A1A2[1]
> A2 = A1A2[2]
>
> FI <- fi(sigma_ADD(estimated_alpha))
> p <- ncol(x)
> EF <- 1 + 1 / p * A2
> DF <- 1 / p^2 * (A1 + 6 * A2 + 2 * p)
> m <- 4 + (p + 2) / (p * DF / (2 * EF^2) - 1)
> lam <- m / (EF * (m - 2))
>
> beta_egls <-
> crossprod(FI,crossprod(x,crossprod(solve(sigma_ADD(estimated_alpha)),y)))
> beta_covmat_old <- FI + 2*crossprod(FI,crossprod(INNER,FI))
>
> dif = beta_egls - beta_tested
>
> F.kr = as.numeric(lam / p *
> crossprod(dif,crossprod(solve(beta_covmat_old),dif)))
> p.v <- pf(F.kr, df1 = p, df2 = m, lower.tail = FALSE, log.p = FALSE)
>
> list <- list(A1 = A1, A2 = A2, EF = EF, DF = DF, df1 = p, df2 = m, F_old
> = F.kr,
> F_new = F.kr, sig_old = p.v, sig_new = p.v, beta_egls = beta_egls,
> beta_covmat_old = beta_covmat_old, beta_covmat_new =
> beta_covmat_old,model = model)
> list }
>
>
> else if (model == "MIX") {
> inner <- function(i,j,alpha) {
> FI <- fi(sigma_MIX(alpha))
> Pi <- pi_MIX(sigma_MIX(alpha),i)
> Pj <- pi_MIX(sigma_MIX(alpha),j)
> Qij <- qij_MIX(sigma_MIX(alpha),i,j)
> Rij <- rij_MIX(sigma_MIX(alpha),i,j)
> W <- expcovmat_MIX(alpha)
> inner = matrix(nrow = ncol(x), ncol = ncol(x))
> inner = W[i,j]*(Qij-crossprod(Pi,crossprod(FI,Pj))-1/4*Rij)
> return(inner) }
> INNER = matrix(0, nrow = ncol(x), ncol = ncol(x))
> for(i in 1:ncol(z)) { for(j in 1:ncol(z)) INNER = INNER +
> inner(i,j,estimated_alpha) }
> a1a2ij <- function(i,j,alpha) {
> FI <- fi(sigma_MIX(alpha))
> Pi <- pi_MIX(sigma_MIX(alpha),i)
> Pj <- pi_MIX(sigma_MIX(alpha),j)
> W <- expcovmat_MIX(alpha)
> a1a2 = vector(mode = "numeric",length = 2)
> a1a2[1] <- W[i,j]*tr(crossprod(FI,Pi))*tr(crossprod(FI,Pj))
> a1a2[2] <- W[i,j]*tr(crossprod(crossprod(Pi,FI),crossprod(FI,Pj)))
> return(a1a2) }
> A1A2 = vector(mode = "numeric", length = 2)
> for(i in 1:ncol(z)) { for(j in 1:ncol(z)) A1A2 = A1A2 +
> a1a2ij(i,j,estimated_alpha) }
> A1 = A1A2[1]
> A2 = A1A2[2]
>
> BIAS_COR_MIX <- function(alpha) {
> SIGMA <- sigma_MIX(alpha)
> FI <- fi(sigma_MIX(alpha))
>
> W <- expcovmat_MIX(alpha)
>
> IN = matrix(0, nrow = nrow(x), ncol = nrow(x))
>
> for(i in 1:ncol(z)) {
> for(j in 1:ncol(z)) IN = IN + W[i,j]*d2sigma.dalphaij_MIX(i,j) } #
> is symmetric
>
> S <- crossprod(solve(SIGMA),crossprod(IN,solve(SIGMA))) # is
> symmetric
> V <- function(t) {
> ancillary1 <- crossprod(S,dsigma.dalphai_MIX(SIGMA,t)) # is
> symmetric
> ancillary2 <- crossprod(x,crossprod(S,x)) # is symmetric
> v1 <- tr(ancillary1)
> v2 <-
> -2*tr(crossprod(sx(SIGMA),crossprod(ancillary1,crossprod(t(x),FI))))
> v3 <-
> -1*tr(crossprod(ancillary2,crossprod(FI,crossprod(pi_MIX(SIGMA,t),FI))))
> v <- v1 + v2 + v3
> return(v) }
>
> bias = matrix(0, nrow = ncol(x), ncol = ncol(x))
> alpha = estimated_alpha
> for(t in 1:ncol(z)) {
> for(s in 1:ncol(z))
> bias = bias +
> (-1/4)*W[t,s]*V(t)*crossprod(FI,crossprod(pi_MIX(SIGMA,s),FI)) }
> return(bias) }
>
> FI <- fi(sigma_MIX(estimated_alpha))
> p <- ncol(x)
> EF <- 1 + 1 / p * A2
> DF <- 1 / p^2 * (A1 + 6 * A2 + 2 * p)
> m <- 4 + (p + 2) / (p * DF / (2 * EF^2) - 1)
> lam <- m / (EF * (m - 2))
>
> beta_egls <-
> crossprod(FI,crossprod(x,crossprod(solve(sigma_MIX(estimated_alpha)),y)))
> beta_covmat_old <- FI + 2*crossprod(FI,crossprod(INNER,FI))
> beta_covmat_new <- beta_covmat_old + BIAS_COR_MIX(estimated_alpha)
>
> dif = beta_egls - beta_tested
>
> F.kr_old = as.numeric(lam / p *
> crossprod(dif,crossprod(solve(beta_covmat_old),dif)))
> p.v_old <- pf(F.kr_old, df1 = p, df2 = m, lower.tail = FALSE, log.p =
> FALSE)
> F.kr_new = as.numeric(lam / p *
> crossprod(dif,crossprod(solve(beta_covmat_new),dif)))
> p.v_new <- pf(F.kr_new, df1 = p, df2 = m, lower.tail = FALSE, log.p =
> FALSE)
>
> list <- list(A1 = A1, A2 = A2, EF = EF, DF = DF, df1 = p, df2 = m, F_old
> = F.kr_old,
> F_new = F.kr_new, sig_old = p.v_old, sig_new = p.v_new, beta_egls =
> beta_egls,
> beta_covmat_old = beta_covmat_old, beta_covmat_new =
> beta_covmat_new,model = model)
> list }
>
> else if (model == "EXP") {
> inner <- function(i,j,alpha) {
> FI <- fi(sigma_EXP(alpha))
> Pi <- pi_EXP(sigma_EXP(alpha),i)
> Pj <- pi_EXP(sigma_EXP(alpha),j)
> Qij <- qij_EXP(sigma_EXP(alpha),i,j)
> Rij <- rij_EXP(sigma_EXP(alpha),i,j)
> W <- expcovmat_EXP(alpha)
> inner = matrix(nrow = ncol(x), ncol = ncol(x))
> inner = W[i,j]*(Qij-crossprod(Pi,crossprod(FI,Pj))-1/4*Rij)
> return(inner) }
> INNER = matrix(0, nrow = ncol(x), ncol = ncol(x))
> for(i in 1:ncol(z)) { for(j in 1:ncol(z)) INNER = INNER +
> inner(i,j,estimated_alpha) }
> a1a2ij <- function(i,j,alpha) {
> FI <- fi(sigma_EXP(alpha))
> Pi <- pi_EXP(sigma_EXP(alpha),i)
> Pj <- pi_EXP(sigma_EXP(alpha),j)
> W <- expcovmat_EXP(alpha)
> a1a2 = vector(mode = "numeric",length = 2)
> a1a2[1] <- W[i,j]*tr(crossprod(FI,Pi))*tr(crossprod(FI,Pj))
> a1a2[2] <- W[i,j]*tr(crossprod(crossprod(Pi,FI),crossprod(FI,Pj)))
> return(a1a2) }
> A1A2 = vector(mode = "numeric", length = 2)
> for(i in 1:ncol(z)) { for(j in 1:ncol(z)) A1A2 = A1A2 +
> a1a2ij(i,j,estimated_alpha) }
> A1 = A1A2[1]
> A2 = A1A2[2]
>
> BIAS_COR_EXP <- function(alpha) {
> SIGMA <- sigma_EXP(alpha)
> FI <- fi(sigma_EXP(alpha))
> W <- expcovmat_EXP(alpha)
>
> IN = matrix(0, nrow = nrow(x), ncol = nrow(x))
>
> for(i in 1:ncol(z)) {
> for(j in 1:ncol(z)) IN = IN + W[i,j]*d2sigma.dalphaij_EXP(SIGMA,i,j) }
> # is symmetric
>
> S <- crossprod(solve(SIGMA),crossprod(IN,solve(SIGMA))) # is
> symmetric
> V <- function(t) {
> ancillary1 <- crossprod(S,dsigma.dalphai_EXP(SIGMA,t)) # is
> symmetric
> ancillary2 <- crossprod(x,crossprod(S,x)) # is symmetric
> v1 <- tr(ancillary1)
> v2 <-
> -2*tr(crossprod(sx(SIGMA),crossprod(ancillary1,crossprod(t(x),FI))))
> v3 <-
> -1*tr(crossprod(ancillary2,crossprod(FI,crossprod(pi_EXP(SIGMA,t),FI))))
> v <- v1 + v2 + v3
> return(v) }
>
> bias = matrix(0, nrow = ncol(x), ncol = ncol(x))
>
> for(t in 1:ncol(z)) {
> for(s in 1:ncol(z))
> bias = bias +
> (-1/4)*W[t,s]*V(t)*crossprod(FI,crossprod(pi_EXP(SIGMA,s),FI)) }
> return(bias) }
>
> FI <- fi(sigma_EXP(estimated_alpha))
> p <- ncol(x)
> EF <- 1 + 1 / p * A2
> DF <- 1 / p^2 * (A1 + 6 * A2 + 2 * p)
> m <- 4 + (p + 2) / (p * DF / (2 * EF^2) - 1)
> lam <- m / (EF * (m - 2))
>
> beta_egls <-
> crossprod(FI,crossprod(x,crossprod(solve(sigma_EXP(estimated_alpha)),y)))
> beta_covmat_old <- FI + 2*crossprod(FI,crossprod(INNER,FI))
> beta_covmat_new <- beta_covmat_old + BIAS_COR_EXP(estimated_alpha)
>
> dif = beta_egls - beta_tested
>
> F.kr_old = as.numeric(lam / p *
> crossprod(dif,crossprod(solve(beta_covmat_old),dif)))
> p.v_old <- pf(F.kr_old, df1 = p, df2 = m, lower.tail = FALSE, log.p =
> FALSE)
> F.kr_new = as.numeric(lam / p *
> crossprod(dif,crossprod(solve(beta_covmat_new),dif)))
> p.v_new <- pf(F.kr_new, df1 = p, df2 = m, lower.tail = FALSE, log.p =
> FALSE)
>
>
> list <- list(A1 = A1, A2 = A2, EF = EF, DF = DF, df1 = p, df2 = m, F_old
> = F.kr_old,
> F_new = F.kr_new, sig_old = p.v_old, sig_new = p.v_new, beta_egls =
> beta_egls,
> beta_covmat_old = beta_covmat_old, beta_covmat_new =
> beta_covmat_new,model = model)
> list }
>
> else return(cat("The wrong argument (model) - must be either >ADD<, or
> >MIX<, or >EXP<.\n")) }
>
>
> kenward_roger_both_print <- function(krprocedure) {
>
> A1 = krprocedure$A1; A2 = krprocedure$A2; EF = krprocedure$EF; DF =
> krprocedure$DF
> p = krprocedure$df1; m = krprocedure$df2; F.kr_old = krprocedure$F_old
> F.kr_new = krprocedure$F_new; p.v_new= krprocedure$sig_new; p.v_old =
> krprocedure$sig_old
> beta_egls = krprocedure$beta_egls; beta_covmat_new =
> krprocedure$beta_covmat_new
> beta_covmat_old = krprocedure$beta_covmat_old; beta_egls =
> krprocedure$beta_egls
>
> if (krprocedure$model == "ADD") {
>
> summary1 <-
> cbind(data.frame("Parameter"=c("Estimates")),BETA=t(as.vector(beta_egls)))
> summary2 <- cbind(BETA=beta_covmat_old)
>
> cat("**Parameter estimates of the additive heteroskedasticity
> model**\n")
>
> cat("-----------------------------------------------------------------------------\n")
> print(summary1,row.names=FALSE)
>
> cat("-----------------------------------------------------------------------------\n")
> cat("**Both original 1997 and new 2009 Kenward-Roger procedure**\n")
> cat(" \n")
> cat("**Kenward-Roger small-sample adjusted estimate of the
> covmat(beta)**\n")
> print(summary2,row.names=FALSE)
> cat(" \n")
> cat("Small-sample adjusted F-statistics =",F.kr_old,"\n")
> cat("df1 (nominator) = ",p,", df2 (denominator) = ",m,", Significance
> =",p.v_old,"\n")
>
> cat("-----------------------------------------------------------------------------\n")
> }
>
> else if (krprocedure$model == "MIX") {
>
> summary1 <-
> cbind(data.frame("Parameter"=c("Estimates")),BETA=t(as.vector(beta_egls)))
> summary2_old <- cbind(BETA=beta_covmat_old)
> summary2_new <- cbind(BETA=beta_covmat_new)
>
> cat("**Parameter estimates of the mixed heteroskedasticity model**\n")
>
> cat("-----------------------------------------------------------------------------\n")
> print(summary1,row.names=FALSE)
>
> cat("-----------------------------------------------------------------------------\n")
> cat("**Original 1997 Kenward-Roger procedure**\n")
> cat(" \n")
> cat("**Kenward-Roger small-sample adjusted estimate of the
> covmat(beta)**\n")
> print(summary2_old,row.names=FALSE)
> cat(" \n")
> cat("Small-sample adjusted F-statistics =",F.kr_old,"\n")
> cat("df1 (nominator) = ",p,", df2 (denominator) = ",m,", Significance
> =",p.v_old,"\n")
>
> cat("-----------------------------------------------------------------------------\n")
> cat("**New 2009 Kenward-Roger procedure**\n")
> cat(" \n")
> cat("**Kenward-Roger small-sample adjusted estimate of the
> covmat(beta)**\n")
> print(summary2_new,row.names=FALSE)
> cat(" \n")
> cat("Small-sample adjusted F-statistics =",F.kr_new,"\n")
> cat("df1 (nominator) = ",p,", df2 (denominator) = ",m,", Significance
> =",p.v_new,"\n")
>
> cat("-----------------------------------------------------------------------------\n")
> }
>
> else {
>
> summary1 <-
> cbind(data.frame("Parameter"=c("Estimates")),BETA=t(as.vector(beta_egls)))
> summary2_old <- cbind(BETA=beta_covmat_old)
> summary2_new <- cbind(BETA=beta_covmat_new)
>
> cat("**Parameter estimates of the multiplicative heteroskedasticity
> model**\n")
>
> cat("-----------------------------------------------------------------------------\n")
> print(summary1,row.names=FALSE)
>
> cat("-----------------------------------------------------------------------------\n")
> cat("**Original 1997 Kenward-Roger procedure**\n")
> cat(" \n")
> cat("**Kenward-Roger small-sample adjusted estimate of the
> covmat(beta)**\n")
> print(summary2_old,row.names=FALSE)
> cat(" \n")
> cat("Small-sample adjusted F-statistics =",F.kr_old,"\n")
> cat("df1 (nominator) = ",p,", df2 (denominator) = ",m,", Significance
> =",p.v_old,"\n")
>
> cat("-----------------------------------------------------------------------------\n")
> cat("**New 2009 Kenward-Roger procedure**\n")
> cat(" \n")
> cat("**Kenward-Roger small-sample adjusted estimate of the
> covmat(beta)**\n")
> print(summary2_new,row.names=FALSE)
> cat(" \n")
> cat("Small-sample adjusted F-statistics =",F.kr_new,"\n")
> cat("df1 (nominator) = ",p,", df2 (denominator) = ",m,", Significance
> =",p.v_new,"\n")
>
> cat("-----------------------------------------------------------------------------\n")
> }}
>
> ######## generovanie vektora ypsilonov
> ##################################################
>
> generate_y <- function(x,z,beta,alpha,model) { res =
> as.matrix(rnorm(nrow(x)))
> if (model == "ADD") {
> sigma_add = diag(as.vector(alpha %*% t(z)), nrow=nrow(z), ncol=nrow(z))
> y = crossprod(t(x),beta) + crossprod(sigma_add^(1/2),res)
> return(y) }
> else if (model == "MIX") {
> sigma_mix = diag((as.vector(alpha %*% t(z)))^2, nrow=nrow(z),
> ncol=nrow(z))
> y = crossprod(t(x),beta) + crossprod(sigma_mix^(1/2),res)
> return(y) }
> else if (model == "EXP") {
> sigma_exp = diag(exp(as.vector(alpha %*% t(z))), nrow=nrow(z),
> ncol=nrow(z))
> y = crossprod(t(x),beta) + crossprod(sigma_exp^(1/2),res)
> return(y) }
> else return(cat("The wrong argument (type)- must be either >ADD<, or
> >MIX<, or >EXP<.\n")) }
>
>
> ######## pomocné funkcie pre simulovanie
> ################################################
>
> loglk <- function(simtype) { if (simtype == "ADD") reml.loglik_ADD
> else if (simtype == "MIX") reml.loglik_MIX
> else reml.loglik_EXP }
>
> gradvec <- function(simtype) { if (simtype == "ADD") gradvec_ADD
> else if (simtype == "MIX") gradvec_MIX
> else gradvec_EXP }
>
> hessmat <- function(simtype) { if (simtype == "ADD") hessmat_ADD
> else if (simtype == "MIX") hessmat_MIX
> else hessmat_EXP }
>
> starter <- function(simtype,x,y) { if (simtype == "ADD")
> abs(coefolsEST_ADD(x,y))
> else if (simtype == "MIX") coefolsEST_MIX(x,y)
> else coefolsEST_EXP(x,y) }
>
> existencecheck <- function(object)
> exists(as.character(substitute(object)))
>
>
>
>
>
>
>
>
>
>
> #######################################################################################
>
> x <- as.matrix(read.table("C:/maddala_70.txt", header = F)) # dané
> pevné x
> z <- as.matrix(read.table("C:/maddala_70.txt", header = F)) # dané
> pevné z
>
>
> nsimsgood=5000; beta=c(0.85,0.90); beta_tested=c(0.85,0.90)
> alpha=c(1, 0.50); realtype="ADD"; simtype="ADD"; prob=0.05
>
> heading = matrix(nrow = 1, ncol = 15 + 2*ncol(z) + 2*ncol(x)); i = 0;
> j = 0;
> name1 <-
>
> c("SimNo","realHET","simHET",rep("realALPHA",ncol(z)),rep("estALPHA",ncol(z)))
> name2 <-
> c("estALPHAlikTYPE",rep("realBETA",ncol(x)),rep("estBETA",ncol(x)))
> name3 <- c("A1","A2","EF","DF","df1","df2")
> name4 <-
> c("Fstatisticsold","Fstatisticsnew","Fquantile","inAreaOld","inAreaNew")
> heading[1,] <- c(name1,name2,name3,name4)
> write.table(heading, file="70addadd__goods.csv", append=T, quote=T,
> sep=";", row.names=F, col.names=F)
> write.table(c("Unsuccessful simulations"), file="70addadd__bads.csv",
> append=T, quote=T, sep=";", row.names=F, col.names=F)
>
>
>
>
> repeat {
> y <<- generate_y(x,z,beta,alpha,realtype); x <<- x; z <<- z
> row = matrix(nrow = 1, ncol = 15 + 2*ncol(z) + 2*ncol(x))
>
> NR <-
>
> try(maxNR(loglk(simtype),gradvec(simtype),hessmat(simtype),start=starter(simtype,x,y)),silent=TRUE)
> if (is.list(NR) == TRUE) kr_nr <-
> kenward_roger_both(simtype,NR$estimate,beta_tested,x,y,z)
>
> if (is.list(NR) == TRUE) NRcode <- NR$code else NRcode <- 99
> if (existencecheck(kr_nr) == TRUE) kr_nrF_old <- kr_nr$F_old else
> kr_nrF_old <- -99
> if (existencecheck(kr_nr) == TRUE) kr_nrF_new <- kr_nr$F_new else
> kr_nrF_new <- -99
> if (existencecheck(kr_nr) == TRUE) qff <-
> qf(prob,df1=kr_nr$df1,df2=kr_nr$df2,lower.tail=TRUE) else qff <- NaN
>
> { if ((NRcode == 0 | NRcode == 1) & (kr_nrF_old >= 0 & kr_nrF_new >= 0 &
> is.nan(qff) == FALSE)) { i = i + 1
> A1 <- kr_nr$A1; A2 <- kr_nr$A2; EF <- kr_nr$EF; DF <- kr_nr$DF
> df1 <- kr_nr$df1; df2 <- kr_nr$df2
> Fold <- kr_nr$F_old; Fnew <- kr_nr$F_new
> f <- qf(prob,df1=kr_nr$df1,df2=kr_nr$df2,lower.tail=TRUE)
>
> rowi1 = c(i, realtype, simtype, alpha, NR$estimate,
> c("Nephton-Raphson"), beta)
> rowi2 = c(t(kr_nr$beta_egls), A1, A2, EF, DF, df1, df2)
> rowi3 = c(Fold, Fnew, f, Fold < f, Fnew < f)
> row[1,] = c(rowi1,rowi2,rowi3)
>
>
> write.table(row,file="70addadd__goods.csv",append=T,quote=T,sep=";",row.names=F,col.names=F)
> rm(y, NR, kr_nr, NRcode, kr_nrF_old, kr_nrF_new, qff) }
>
> else { BFGS <-
>
> try(maxBFGS(loglk(simtype),gradvec(simtype),hessmat(simtype),start=starter(simtype,x,y)),silent=TRUE)
> if (is.list(BFGS) == TRUE) kr_bfgs <-
> kenward_roger_both(simtype,BFGS$estimate,beta_tested,x,y,z)
>
> if (is.list(BFGS) == TRUE) BFGScode <- BFGS$code else BFGScode <- 99
> if (existencecheck(kr_bfgs) == TRUE) kr_bfgsF_old <- kr_bfgs$F_old else
> kr_bfgsF_old <- -99
> if (existencecheck(kr_bfgs) == TRUE) kr_bfgsF_new <- kr_bfgs$F_new else
> kr_bfgsF_new <- -99
> if (existencecheck(kr_bfgs) == TRUE) qff <-
> qf(prob,df1=kr_bfgs$df1,df2=kr_bfgs$df2,lower.tail=TRUE) else qff <- NaN
>
> { if ((BFGScode == 0) & (kr_bfgsF_old >= 0 & kr_bfgsF_new >= 0 &
> is.nan(qff) == FALSE)) { i = i + 1
> A1 <- kr_bfgs$A1; A2 <- kr_bfgs$A2; EF <- kr_bfgs$EF; DF <- kr_bfgs$DF
> df1 <- kr_bfgs$df1; df2 <- kr_bfgs$df2
> Fold <- kr_bfgs$F_old; Fnew <- kr_bfgs$F_new
> f <- qf(prob,df1=kr_bfgs$df1,df2=kr_bfgs$df2,lower.tail=TRUE)
>
> rowi1 = c(i, realtype, simtype, alpha, BFGS$estimate, c("BFGS"), beta)
> rowi2 = c(t(kr_bfgs$beta_egls), A1, A2, EF, DF, df1, df2)
> rowi3 = c(Fold, Fnew, f, Fold < f, Fnew < f)
> row[1,] = c(rowi1,rowi2,rowi3)
>
>
> write.table(row,file="70addadd__goods.csv",append=T,quote=T,sep=";",row.names=F,col.names=F)
> rm(y, BFGS, kr_bfgs, BFGScode, kr_bfgsF_old, kr_bfgsF_new, qff) }
>
> else { NM <-
>
> try(maxNM(loglk(simtype),gradvec(simtype),hessmat(simtype),start=starter(simtype,x,y)),silent=TRUE)
> if (is.list(NM) == TRUE) kr_nm <-
> kenward_roger_both(simtype,NM$estimate,beta_tested,x,y,z)
>
> if (is.list(NM) == TRUE) NMcode <- NM$code else NMcode <- 99
> if (existencecheck(kr_nm) == TRUE) kr_nmF_old <- kr_nm$F_old else
> kr_nmF_old <- -99
> if (existencecheck(kr_nm) == TRUE) kr_nmF_new <- kr_nm$F_new else
> kr_nmF_new <- -99
> if (existencecheck(kr_nm) == TRUE) qff <-
> qf(prob,df1=kr_nm$df1,df2=kr_nm$df2,lower.tail=TRUE) else qff <- NaN
>
> { if ((NMcode == 0) & (kr_nmF_old >= 0 & kr_nmF_new >= 0 & is.nan(qff)
> == FALSE)) { i = i + 1
> A1 <- kr_nm$A1; A2 <- kr_nm$A2; EF <- kr_nm$EF; DF <- kr_nm$DF
> df1 <- kr_nm$df1; df2 <- kr_nm$df2
> Fold <- kr_nm$F_old; Fnew <- kr_nm$F_new
> f <- qf(prob,df1=kr_nm$df1,df2=kr_nm$df2,lower.tail=TRUE)
>
> rowi1 = c(i, realtype, simtype, alpha, NM$estimate, c("Nelder-Mead"),
> beta)
> rowi2 = c(t(kr_nm$beta_egls), A1, A2, EF, DF, df1, df2)
> rowi3 = c(Fold, Fnew, f, Fold < f, Fnew < f)
> row[1,] = c(rowi1,rowi2,rowi3)
>
>
> write.table(row,file="70addadd__goods.csv",append=T,quote=T,sep=";",row.names=F,col.names=F)
> rm(y, NN, kr_nm, NMcode, kr_nmF_old, kr_nmF_new, qff) }
>
> else { SANN <-
>
> try(maxSANN(loglk(simtype),gradvec(simtype),hessmat(simtype),start=starter(simtype,x,y)),silent=TRUE)
> if (is.list(SANN) == TRUE) kr_sann <-
> kenward_roger_both(simtype,SANN$estimate,beta_tested,x,y,z)
>
> if (is.list(SANN) == TRUE) SANNcode <- SANN$code else SANNcode <- 99
> if (existencecheck(kr_sann) == TRUE) kr_sannF_old <- kr_sann$F_old else
> kr_sannF_old <- -99
> if (existencecheck(kr_sann) == TRUE) kr_sannF_new <- kr_sann$F_new else
> kr_sannF_new <- -99
> if (existencecheck(kr_sann) == TRUE) qff <-
> qf(prob,df1=kr_sann$df1,df2=kr_sann$df2,lower.tail=TRUE) else qff <- NaN
>
> { if ((SANNcode == 0) & (kr_sannF_old >= 0 & kr_sannF_new >= 0 &
> is.nan(qff) == FALSE)) { i = i + 1
> A1 <- kr_sann$A1; A2 <- kr_sann$A2; EF <- kr_sann$EF; DF <- kr_sann$DF
> df1 <- kr_sann$df1; df2 <- kr_sann$df2
> Fold <- kr_sann$F_old; Fnew <- kr_sann$F_new
> f <- qf(prob,df1=kr_sann$df1,df2=kr_sann$df2,lower.tail=TRUE)
>
> rowi1 = c(i, realtype, simtype, alpha, SANN$estimate, c("SANN"), beta)
> rowi2 = c(t(kr_sann$beta_egls), A1, A2, EF, DF, df1, df2)
> rowi3 = c(Fold, Fnew, f, Fold < f, Fnew < f)
> row[1,] = c(rowi1,rowi2,rowi3)
>
>
> write.table(row,file="70addadd__goods.csv",append=T,quote=T,sep=";",row.names=F,col.names=F)
> rm(y, SANN, kr_sann, SANNcode, kr_sannF_old, kr_sannF_new, qff) }
>
> else { j = j + 1
>
> write.table(j,file="70addadd__bads.csv",append=T,quote=T,sep=";",row.names=F,col.names=F)
> rm(y, NR, BFGS, NM, SANN, NRcode, BFGScode, NMcode, SANNcode,
> kr_nrF_old, kr_bfgsF_old, kr_nmF_old, kr_sannF_old)
> rm(qff, kr_nrF_new, kr_bfgsF_new, kr_nmF_new, kr_sannF_new)
> if (existencecheck(kr_nr) == TRUE) rm(kr_nr)
> if (existencecheck(kr_bfgs) == TRUE) rm(kr_bfgs)
> if (existencecheck(kr_nm) == TRUE) rm(kr_nm)
> if (existencecheck(kr_sann) == TRUE) rm(kr_sann) } } } } } } } } }
>
> if (i == nsimsgood) break
>
>
>
> cat("-----------------------------------------------------------------------------\n")
> cat("Number of successful good simulations = ",i,"\n")
> cat("Number of unsuccessful simulations = ",j,"\n")
>
> cat("-----------------------------------------------------------------------------\n")
>
> ______________________________________________
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html<http://www.r-project.org/posting-guide.html>
> and provide commented, minimal, self-contained, reproducible code.
>
--
Jim Holtman
Cincinnati, OH
+1 513 646 9390
What is the problem that you are trying to solve?
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
