################################# #### 2022-11-18, Friday ##### ################################# ################################# #### 2022-12-19, Monday ##### ################################# # I debugged a bug in f.ll, d.ll.theta, d2.ll.theta, Prod1Derivative ################################## #### 2023-01-10, Tuesday ##### ################################## # I added a new function Level2.EM.0 #################################### #### 2023-01-11, Wednesday ##### #################################### # I added a new function SynMP.0 #################################### #### 2023-03-22,23 ##### #################################### # I added a new input argument to the function SynMP.0. Users can now increase the iteration of the EM algorithm to have more accurate starting # values for the Newton method. # Remove a bug in SynMP.0: the bug is related to SE estimates #################################### #### 2023-08-03, Thursday ##### #################################### # I added a new input argument to the function SynMP.0. Users can now specify the estimation to be the EM algorithm and it will compute SEs accordingly. # I added two new functions Make.Select.Matrix() and Fixed.wls(). These two functions faciliate the difference test of a reduced model against the full model. #################################### #### 2024-05-02, Thursday ##### #################################### # I added a new function Select.A to select AR weights from the first stage analysis. # The file contains ten R functions. # f.ll <- function(Beta, D, Beta.star.i, CMatrix.i,TotalA) computes the loglikelihood function # d.ll.theta <- function(Beta, D, Beta.star.i, CMatrix.i, TotalA) computes the derivatives of the loglikelihood WRT parameters # Level2.EM.0 <- function (BetaStar,TotalDesign, TotalCMatrix, nIteration=15) performs the EM-like algorithm # Level2.Newton.0 <- function (theta.0, BetaStar,TotalDesign, TotalCMatrix, nIteration=15, Newton.epsilon=0.0001) estimates the model using the Newton-Raphson method # d2.ll.theta <- function (BetaStar,TotalDesign, TotalCMatrix, Fixed, Cov.Random) computes the second order derivatives of LLF WRT parameters # Prod1Derivative <- function(Beta, D, Beta.star.i, CMatrix.i, TotalA) computes the product of the first order derivatives # SynMP.0 <- function(BetaStar,TotalDesign,TotalCMatrix) synthesizes time series of multiple persons # Make.Select.Matrix <- function(log.select) makes a selection matrix according to a vector of logical variables # Fixed.wls <- function(BetaStar,TotalDesign,TotalCMatrix) performs wls on variables of fixed effects # Select.A <- function (filename) selects AR weights from the first stage analysis. #---------------------------------------------------------------------------------- f.ll <- function(Beta, D, Beta.star.i, CMatrix.i,TotalA) { n = dim(Beta.star.i)[1] p = dim(Beta.star.i)[2] loglikelihood = 0 for (i in 1:n) { term1 = log(det(CMatrix.i[,,i]+D)) term2 = sum( solve((CMatrix.i[,,i]+D),(Beta.star.i[i,] - TotalA[,,i]%*% Beta)) * (Beta.star.i[i,] - TotalA[,,i]%*% Beta)) loglikelihood = loglikelihood + term1 + term2 } # return(loglikelihood) return(loglikelihood*(-0.5)) # A bug found on 2022-12-19 } # function f.ll #-------------------------------------------------------------------------------------- d.ll.theta <- function(Beta, D, Beta.star.i, CMatrix.i, TotalA) { itemp = dim(TotalA) p = itemp[1] m = itemp[2] n = itemp[3] d.ll.2.Beta <- rep(0,m) for (i in 1:n) { d.ll.2.Beta = d.ll.2.Beta + t(TotalA[,,i]) %*% solve((CMatrix.i[,,i]+D),(Beta.star.i[i,] - TotalA[,,i]%*% Beta)) } d.ll.2.D <- matrix(0,p,p) for (i in 1:n) { Temp.M.1 = solve(CMatrix.i[,,i]+D) Temp.v.1 = array( solve((CMatrix.i[,,i]+D),(Beta.star.i[i,] - TotalA[,,i]%*% Beta))) # Make it to an array d.ll.2.D = d.ll.2.D + Temp.M.1 - Temp.v.1 %o% Temp.v.1 } # d.ll.2.D Temp.Sym = d.ll.2.D + t(d.ll.2.D) - diag(diag(d.ll.2.D)) # return(list(d.f.Beta = d.ll.2.Beta * (-2), d.f.D = d.ll.2.D, d.f.theta = c(d.ll.2.Beta * (-2), # array(Temp.Sym[upper.tri(Temp.Sym, diag=TRUE)])) )) return(list(d.f.Beta = d.ll.2.Beta, d.f.D = d.ll.2.D * (-0.5), d.f.theta = c(d.ll.2.Beta, (-0.5) * array(Temp.Sym[upper.tri(Temp.Sym, diag=TRUE)])) )) # A bug found on 2022-12-19 } # function d.ll.theta #---------------------------------------------------------------------------- # The function Level2.EM.0 synthesizes parameter estimates from different individuals. # The algorithm is a EM-like method. # 2022-11-17, Thursday! Guangjian Zhang Level2.EM.0 <- function (BetaStar,TotalDesign, TotalCMatrix, nIteration=15) { # The function has four input arguments. # BetaStar, n.sample:p, the matrix of individual parameter estimates # TotalDesign, p:m:n.sample, the design matrix of each individual # TotalCMatrix, p:p:n.sample, the estimation error covariance matrix of each individual # nIteration, a scalar, the maximum number of iterations ### Step 0, housekeeping Beta.star.i = BetaStar TotalA = TotalDesign CMatrix.i = TotalCMatrix n.iteration = nIteration itemp = dim(TotalA) p = itemp[1] m = itemp[2] n.sample = itemp[3] ### Step 1, starting values term1 = rep(0,m) term2 = matrix(0,m,m) for (i in 1:n.sample) { # term1 = term1 + t(TotalA[,,i]) %*% Temp[i,] term1 = term1 + t(TotalA[,,i]) %*% Beta.star.i[i,] # Remove a bug 2022-12-21, GZ term2 = term2 + t(TotalA[,,i]) %*% TotalA[,,i] } Beta.current = solve(term2,term1) residual = matrix(NA,n.sample,p) for (i in 1:n.sample) { # residual[i,] = Temp[i,] - TotalA[,,i] %*% Beta.current residual[i,] = Beta.star.i[i,] - TotalA[,,i] %*% Beta.current # Remove a bug 2022-12-21, GZ } D.current = t(residual) %*% residual / n.sample ### Step 2, the EM-like algorithm Beta.iteration = matrix(NA, n.iteration, m) D.iteration = array(rep(NA, p*p*n.iteration), dim=c(p,p,n.iteration)) Beta.i.current = matrix(NA, n.sample,p) Term1 = matrix(NA,n.sample,p) # I made a mistake at the beginning: it should be p for (i in 1:n.sample) { Term1[i,] = solve(CMatrix.i[,,i],Beta.star.i[i,]) } for (k in 1:n.iteration) { # The E step for (i in 1:n.sample) { Test = solve(D.current, TotalA[,,i]%*% Beta.current) Beta.i.current [i,] = solve( (solve(CMatrix.i[,,i]) + solve(D.current)) , (Term1[i,] + Test) ) } # The M Step Temp = matrix(0,m,m) for (i in 1:n.sample) { Temp = Temp + t(TotalA[,,i]) %*% solve(D.current, TotalA[,,i]) # The relationship between D and A is strange! } Temp.v = rep(0,m) for (i in 1:n.sample) { Temp.v = Temp.v + t(TotalA[,,i]) %*% solve(D.current, Beta.i.current [i,]) } Beta.next = solve(Temp, Temp.v) ## update D Temp = matrix(0,p,p) for (i in 1:n.sample) { Temp.v = Beta.i.current [i,] - TotalA[,,i] %*% Beta.next Temp = Temp + solve(solve(CMatrix.i[,,i]) + solve(D.current)) + Temp.v %*% t(Temp.v) } D.next = Temp / n.sample Beta.iteration[k,] = Beta.next D.iteration[,,k] = D.next # Update Beta and D Beta.current = Beta.next D.current = D.next } # (k in 1:n.iteration) ### output the results return(list(BetaAll=Beta.iteration, DAll = D.iteration)) } # Level2.EM.0 #---------------------------------------------------------------------------------------------------- # The function Level2.Newton.0 synthesizes parameter estimates from different individuals. # The algorithm is the Newton-Raphson method # 2023-01-10, Tuesday! Guangjian Zhang Level2.Newton.0 <- function (theta.0, BetaStar,TotalDesign, TotalCMatrix, nIteration=15, Newton.epsilon=0.0001) { # arguments # The function has five input arguments. # theta.0, m+p(p+1)/2, starting values for the parameters # BetaStar, n.sample:p, the matrix of individual parameter estimates # TotalDesign, p:m:n.sample, the design matrix of each individual # TotalCMatrix, p:p:n.sample, the estimation error covariance matrix of each individual # nIteration, a scalar, the maximum number of iterations # Newton.epsilon, a scalar, the tolerance of the largest first order derivative, convergence criterion # Step 0, housekeeping itemp = dim(TotalDesign) p = itemp[1] m = itemp[2] n.sample = itemp[3] convergence = FALSE # Step 1, prepare starting values theta.t <- theta.0 beta.t <- theta.t[1:m] D.t <- matrix(NA,p,p) vecD <- theta.t[(m+1):(m+p*(p+1)/2)] m1 <- matrix(NA,p,p) m1[upper.tri(m1, diag=TRUE)] <- vecD D.t = t(m1) D.t[upper.tri(D.t, diag=TRUE)] <- vecD # step 2, iterations for (t in 1:nIteration) { # Check convergence Derivative1.t = d.ll.theta(Beta = beta.t, D=D.t, Beta.star.i=BetaStar, CMatrix.i=TotalCMatrix, TotalA=TotalDesign) if (max(abs(Derivative1.t$d.f.theta))j) { ijs = i*(i-1)/2 + j } else { ijs = j*(j-1)/2 + i } M.Select[ij,ijs] = 1 } # i } # j D2theta[1:m,1:m] = d2.Beta.Beta D2theta[1:m, (m + 1) : (m + p*(p+1)/2) ] = d2.Beta.D D2theta[(m + 1) : (m + p*(p+1)/2) , 1:m ] = t(d2.Beta.D) D2theta[(m + 1) : (m + p*(p+1)/2), (m + 1) : (m + p*(p+1)/2)] = t(M.Select) %*% d2.D.D %*% M.Select # return(list(Full=D2theta, BetaBeta=d2.Beta.Beta, BetaD=d2.Beta.D, DD=d2.D.D)) # A bug found on 2022-12-09 return(list(Full=D2theta * (-0.5), BetaBeta=d2.Beta.Beta * (-0.5), BetaD=d2.Beta.D * (-0.5), DD=d2.D.D * (-0.5))) } # d2.ll.theta #----------------------------------------------------------------------------------------------------- Prod1Derivative <- function(Beta, D, Beta.star.i, CMatrix.i, TotalA) { itemp = dim(TotalA) p = itemp[1] m = itemp[2] n = itemp[3] Derivative.Beta <- matrix(0,n,m) for (i in 1:n) { Derivative.Beta[i, 1:m] = t(TotalA[,,i]) %*% solve((CMatrix.i[,,i]+D),(Beta.star.i[i,] - TotalA[,,i]%*% Beta)) } Derivative.Beta = Derivative.Beta * (-2) Derivative.D <- matrix(0,n, p*(p+1)/2) for (i in 1:n) { Temp.M.1 = solve(CMatrix.i[,,i]+D) Temp.v.1 = array( solve((CMatrix.i[,,i]+D),(Beta.star.i[i,] - TotalA[,,i]%*% Beta))) # Make it to an array Temp2Dim = Temp.M.1 - Temp.v.1 %o% Temp.v.1 Temp.Sym = Temp2Dim + t(Temp2Dim) - diag(diag(Temp2Dim)) Derivative.D[i, 1:(p*(p+1)/2)] = array(Temp.Sym[upper.tri(Temp.Sym, diag=TRUE)]) } Product = matrix(0, (m + p*(p+1)/2), (m + p*(p+1)/2)) Product[1:m, 1:m] = t(Derivative.Beta) %*% Derivative.Beta Product[1:m, (m+1):(m+p*(p+1)/2)] = t(Derivative.Beta) %*% Derivative.D Product[(m+1):(m+p*(p+1)/2), 1:m] = t(Derivative.D) %*% Derivative.Beta Product[(m+1):(m+p*(p+1)/2), (m+1):(m+p*(p+1)/2)] = t(Derivative.D) %*% Derivative.D # return(list(Product = Product, Derivative.Beta = Derivative.Beta, Derivative.D = Derivative.D )) # A bug found on 2022-12-19, Guangjian return(list(Product = Product/4, Derivative.Beta = Derivative.Beta/4, Derivative.D = Derivative.D/4 )) } # function Prod1Derivative #------------------------------------------------------------------------------------------------------ # The function SynMP.0 wraps several functions for synthesize time series from multiple persons. # 2023-03-22, Thursday! Guangjian Zhang. Adding the input argument of estimation = EM. # 2023-03-22, Wednesday! Guangjian Zhang. Adding the facility of Iter.EM to have more accurate starting values # 2023-01-11, Wednesday! Guangjian Zhang. SynMP.0 <- function(BetaStar,TotalDesign,TotalCMatrix, Iter.EM = 10, estimation = 'NewtonRaphson') { # arguments # The function has three input arguments. # BetaStar, n.sample:p, the matrix of individual parameter estimates # TotalDesign, p:m:n.sample, the design matrix of each individual # TotalCMatrix, p:p:n.sample, the estimation error covariance matrix of each individual # It directly invokes four other functions: Level2.EM.0, Level2.Newton.0, Prod1Derivative, and d2.ll.theta # step 0, housekeeping itemp = dim(TotalDesign) p = itemp[1] m = itemp[2] n.sample = itemp[3] # step 1, point estimates Starting = Level2.EM.0(BetaStar=BetaStar,TotalDesign=TotalDesign, TotalCMatrix=TotalCMatrix, nIteration=Iter.EM) beta.t = Starting$BetaAll[Iter.EM,] D.t = Starting$DAll[,,Iter.EM] theta.t <- c(beta.t,D.t[upper.tri(D.t, diag=TRUE)]) convergence = FALSE I.convergence = -1 if (estimation == 'NewtonRaphson') { Test.Newton = Level2.Newton.0 (theta.0=theta.t, BetaStar=BetaStar,TotalDesign=TotalDesign, TotalCMatrix=TotalCMatrix) beta.Newton = Test.Newton$theta[1:m] vecD <- Test.Newton$theta[(m+1):(m+p*(p+1)/2)] m1 <- matrix(NA,p,p) m1[upper.tri(m1, diag=TRUE)] <- vecD D.Newton = t(m1) D.Newton[upper.tri(D.Newton, diag=TRUE)] <- vecD theta.Newton = Test.Newton$theta convergence = Test.Newton$convergence I.convergence = Test.Newton$conI } else { beta.Newton = beta.t D.Newton = D.t theta.Newton = c(beta.t,D.t[upper.tri(D.t, diag=TRUE)]) } # step 2, standard errors TestProduct = Prod1Derivative(Beta=beta.Newton, D=D.Newton, Beta.star.i=BetaStar, CMatrix.i=TotalCMatrix, TotalA=TotalDesign) Test2Derivative = d2.ll.theta(BetaStar=BetaStar,TotalDesign=TotalDesign, TotalCMatrix=TotalCMatrix, Fixed = beta.Newton, Cov.Random=D.Newton) acov.info = solve(-Test2Derivative$Full) # Second order derivatives, remember to add the negative sign acov.sand = solve(Test2Derivative$Full) %*% (TestProduct$Product) %*% solve(Test2Derivative$Full) # The sandwich SE estimator # SE.info = sqrt(diag(acov.info) / n.sample ) # It is wrong to rescale with n.sample, because the parameter variance estimate # SE.sand = sqrt(diag(acov.sand) / n.sample ) # has already been scaled. All the terms are the sum of n.sample terms. # It is different from covariance structure modeling with the covariance matrices SE.info = sqrt( diag(acov.info) ) # or correlation matrices. As n.sample increases, these matrices converge to the population SE.sand = sqrt( diag(acov.sand) ) # counterpart. The bug was found and removed on 2023-03-023, Thursday, Guangjian Zhang. # outputs: point estimate, standard errors, the loglikelihood function, and more outputs LogL = f.ll(Beta=beta.Newton, D=D.Newton, Beta.star.i=BetaStar, CMatrix.i=TotalCMatrix,TotalA=TotalDesign) CI.level = 0.95 CI.L.info <- theta.Newton + qnorm((1-CI.level)/2) * SE.info CI.U.info <- theta.Newton + qnorm((1+CI.level)/2) * SE.info CI.L.sand <- theta.Newton + qnorm((1-CI.level)/2) * SE.sand CI.U.sand <- theta.Newton + qnorm((1+CI.level)/2) * SE.sand SE.beta.sand = SE.sand[1:m] vecD <- SE.sand[(m+1):(m+p*(p+1)/2)] m1 <- matrix(NA,p,p) m1[upper.tri(m1, diag=TRUE)] <- vecD SE.D.sand = t(m1) SE.D.sand[upper.tri(D.Newton, diag=TRUE)] <- vecD more = list(convergence = convergence, I.convergence = I.convergence, SE.info=SE.info, SE.sand=SE.sand, CI.L.info=CI.L.info, CI.U.info=CI.U.info,CI.L.sand=CI.L.sand,CI.U.sand=CI.U.sand) return(list(beta = beta.Newton, Sigma = D.Newton, SE.beta = SE.beta.sand, SE.Sigma = SE.D.sand, LogL = LogL, more = more, estimation = estimation )) } # SynMP.0 #---------------------------------------------------------------------------------- Make.Select.Matrix <- function(log.select) { ## make a selection matrix according to the vector of logical variables ## An element is selected if the corresponding log variable is TRUE i.full = length(log.select) i.select = which(log.select) i.small = length(i.select) M.select <- matrix(0,i.small,i.full) for (i in 1:i.small) { M.select[i,i.select[i]] = 1 } return(M.select) } # Make.Select.Matrix #--------------------------------------------------------------------------------- Fixed.wls <- function(BetaStar,TotalDesign,TotalCMatrix) { ## Perform weighted least squares estimation for the intercepts and the slopes ## of the fixed effects of the inter-invidual stage itemp = dim(TotalDesign) p = itemp[1] m = itemp[2] n = itemp[3] Term1 = matrix(0,m,m) Term2 = matrix(0,m,1) for (i in 1:n) { Term1 = Term1 + t(TotalDesign[,,i]) %*% solve(TotalCMatrix[,,i], TotalDesign[,,i]) Term2 = Term2 + t(TotalDesign[,,i]) %*% solve(TotalCMatrix[,,i], BetaStar[i,1:p]) } beta = solve(Term1, Term2) SE = sqrt(diag(solve(Term1))) return(list(beta=beta, SE = SE)) } # Fixed.wls #------------------------------------------------------------------------------------ Select.A <- function (filename) { # The function Select.A selects the point estimates and their estimation error matrix # from two files produced by DyFA1L # One needs to change the values mv,nf, nar, nma, nlag, isample, and bnumber according to different # Model specifications # Assign parameter values according to your model mv = 7 nf=2 nar = 1 nma = 0 nlag = 1 isample = 81 bnumber = 200 ntotal = 10 + mv*nf + mv + nf*nf + nar*nf*nf + nma*nf*nf + (nlag + 1)*nf*nf + nf*nf # Read in sample data and bootstrap data Sample.point = array(scan(paste0(trimws(filename),'.spk')), dim=c(ntotal,isample)) Boot.point = array(scan(paste0(trimws(filename),'.bpk')), dim=c(ntotal,isample*bnumber)) # Determine the order by which the parameters are presented Lambda.l = 11 Lambda.u = 10 + mv*nf epsilon.l = Lambda.u + 1 epsilon.u = Lambda.u + mv Psi.l = epsilon.u + 1 Psi.u = epsilon.u + nf*nf A.l = Psi.u + 1 A.u = Psi.u + nar*nf*nf B.l = A.u + 1 B.u = A.u + nma*nf*nf Phi.l = Psi.u + (nar+nma)*nf*nf + 1 Phi.u = Psi.u + (nar+nma)*nf*nf + (nlag+1)*nf*nf Theta.l = Phi.u + 1 Theta.u = Phi.u + nf*nf # Select parameter values of AR weights Sample.A = Sample.point[A.l:A.u,1:isample] Boot.A = Boot.point[A.l:A.u,1:(isample*bnumber)] # Compute estimation error VError = array(rep(NA,nar*nf*nf*nar*nf*nf*isample),dim=c(nar*nf*nf,nar*nf*nf,isample)) for (i in 1:isample) { itemp = (i - 1) * bnumber + 1 VError[1:(nar*nf*nf), 1:(nar*nf*nf), i] = cov(t(Boot.A[1:(nar*nf*nf),itemp:(i*bnumber)])) } # Returning the point estimates and estimation error return(list(Point = t(Sample.A), VError=VError)) } # Select.A <- function (filename) #---------------------------------------------------------------------------------------------