// ASSUMES Rparameterization CODES // 1 alpha, tau // 2 alpha, omega // 3 eta, tau // 4 eta, omega // 5 iota, tau // 6 iota, omega // 7 eta, kappa // 8 phi // 9 ln phi capture mata: mata drop xtmixediou_estimation() capture mata: mata drop NR_REML_d2() capture mata: mata drop crossProductMatrices() capture mata: mata drop sweep_partitionK() capture mata: mata drop Gdot_logcholG() capture mata: mata drop Gdot_arctanG() capture mata: mata drop Gdotdot_logcholG() capture mata: mata drop Gdotdot_arctanG() capture mata: mata drop IOUcov() capture mata: mata drop Browniancov() capture mata: mata drop Rdot_at() capture mata: mata drop Rdotdot_at() capture mata: mata drop Rdot_ao() capture mata: mata drop Rdotdot_ao() capture mata: mata drop Rdot_et() capture mata: mata drop Rdotdot_et() capture mata: mata drop Rdot_eo() capture mata: mata drop Rdotdot_eo() capture mata: mata drop Rdot_it() capture mata: mata drop Rdotdot_it() capture mata: mata drop Rdot_io() capture mata: mata drop Rdotdot_io() capture mata: mata drop Rdot_ek() capture mata: mata drop Rdotdot_ek() capture mata: mata drop Rdot_phi() capture mata: mata drop Rdotdot_phi() capture mata: mata drop Rdot_lnphi() capture mata: mata drop Rdotdot_lnphi() capture mata: mata drop Rconversion() capture mata: mata drop tabulatingBeta() capture mata: mata drop tabTheta_sigmaSq() mata: void xtmixediou_estimation(string vector id, string vector dependent, string vector fevars, string vector reffects, string vector time, /// real scalar Rparameterization, string rowvector initialValues, string matrix scheme, string scalar singularHmethod, /// string scalar tr_log, string scalar tr_trace, string scalar tr_gradient, string scalar tr_showstep, string scalar tr_hessian, /// real scalar maxIterations) { real scalar todo, lnf, sigmaSquared, converged, numberOfIterations, colsZ, numberOfGparas, ll_reml, iter_algorithm, type_algorithm real scalar errorCode, returnCode, IOU1, IOU2, numberOfPanels, numberOfObs, min_ni, max_ni, numFE, cols_b, isPositiveDefinite real scalar iter0start, final_sigma, final_sigmaSquared, final_alpha, final_tau, final_phi string scalar errorText real rowvector startingValues, optimizedParameters, settings, gradient, b, final_Rparameters real colvector Y, idVar, timeVar, beta, final_beta real matrix Hessian, X, Z, panelInfo, varbeta, final_varbeta, table_beta, table_theta, V, algorithm real matrix final_G, final_arctanG transmorphic S // CREATE VIEWS st_view(idVar=., ., id) st_view(Y=., ., dependent) st_view(X=., ., tokens(fevars)) st_view(Z=., ., tokens(reffects)) st_view(timeVar=., ., tokens(time)) startingValues = st_matrix(initialValues) algorithm = st_matrix(scheme) // CONVERT THE startingValues TO THE CORRECT PARAMETERIZATION BEFORE PERFORMING OPTIMIZATION // IF THE IOU MODEL THE startingValues WILL ALWAYS CONTAIN logcholGstar, aLpha AND taustar // IF THE BROWNIAN MOTION MODEL THE startingValues WILL ALWAYS CONTAIN logcholGstar, phistar colsZ = cols(Z) numberOfGparas = 0.5*colsZ*(colsZ+1) optimizedParameters = startingValues if (Rparameterization > 1 & Rparameterization < 8) { // IOU MODEL Rconversion(1, startingValues[1,numberOfGparas+1], startingValues[1,numberOfGparas+2], /// Rparameterization, IOU1=., IOU2=.) optimizedParameters[1,numberOfGparas+1] = IOU1 optimizedParameters[1,numberOfGparas+2] = IOU2 } // PANEL SET-UP panelInfo = panelsetup(idVar,1) numberOfPanels = panelstats(panelInfo)[1] numberOfObs = panelstats(panelInfo)[2] min_ni = panelstats(panelInfo)[3] max_ni = panelstats(panelInfo)[4] // SETS ALL PRINTED RESULTS TO THE SCREEN AS TEXT displayas("text") /******************************** ESTIMATION ********************************/ printf("Performing Hessian-based optimization\n") iter_algorithm = 1 converged = 0 iter0start = 1 // TO ENABLE TO ENTER WHILE LOOP WHEN iterate(0) IS SPECIFIED numberOfIterations = 0 while(converged == 0 & (iter0start == 1 | numberOfIterations < maxIterations)) { type_algorithm = algorithm[iter_algorithm,1] settings = (1,Rparameterization,type_algorithm) maxIter_algorithm = algorithm[iter_algorithm,2] // IF THE NEXT maxIter_algorithm WILL EXCEED THE TOTAL NUMBER OF ITERATIONS THEN RESET maxIter_algorithm if (numberOfIterations+maxIter_algorithm >= maxIterations) { maxIter_algorithm = maxIterations - numberOfIterations - 1 } if (rows(algorithm) > 1) { if (type_algorithm ==1) { if (iter0start==1) { printf("(setting algorithm to nr)\n") } else { printf("(switching algorithm to nr)\n") } } else if (type_algorithm ==2) { if (iter0start==1) { printf("(setting algorithm to fs)\n") } else { printf("(switching algorithm to fs)\n") } } else { if (iter0start==1) { printf("(setting algorithm to ai)\n") } else { printf("(switching algorithm to ai)\n") } } } iter0start = 0 // SWITCH OFF THEREAFTER S = optimize_init() optimize_init_which(S, "min") optimize_init_evaluator(S, &NR_REML_d2()) optimize_init_conv_maxiter(S, maxIter_algorithm) optimize_init_evaluatortype(S, "d2") optimize_init_params(S, optimizedParameters) optimize_init_argument(S, 1, X) optimize_init_argument(S, 2, Y) optimize_init_argument(S, 3, Z) optimize_init_argument(S, 4, panelInfo) optimize_init_argument(S, 5, timeVar) optimize_init_argument(S, 6, settings) optimize_init_argument(S, 7, sigmaSquared) optimize_init_argument(S, 8, beta) optimize_init_argument(S, 9, varbeta) optimize_init_valueid(S, "-2 x log restricted-likelihood") optimize_init_conv_warning(S, "off") optimize_init_singularHmethod(S, singularHmethod) optimize_init_tracelevel(S, tr_log) optimize_init_trace_params(S, tr_trace) optimize_init_trace_gradient(S, tr_gradient) optimize_init_trace_Hessian(S, tr_hessian) optimize_init_trace_step(S,tr_showstep) errorCode = _optimize(S) if (errorCode > 0) { returnCode = optimize_result_returncode(S) errorText = optimize_result_errortext(S) printf(errorText) exit(error(returnCode)) } converged = optimize_result_converged(S) numberOfIterations = numberOfIterations + optimize_result_iterations(S) optimizedParameters = optimize_result_params(S) iter_algorithm = iter_algorithm + 1 if(iter_algorithm > rows(algorithm)) iter_algorithm = 1 } // END OF while LOOP ll_reml = optimize_result_value(S) ll_reml = -0.5*ll_reml final_sigmaSquared = sigmaSquared final_sigma = sqrt(final_sigmaSquared) if(converged == 0) { printf("Convergence not achieved\n") if(numberOfIterations >= maxIterations) { errorCode = 3360 } } /***************************************** CALCULATION OF STANDARD ERRORS *****************************************/ // CONVERT PARAMETERS FROM logcholGstar -> G final_G = invvech(optimizedParameters[|1,1 \ 1,numberOfGparas|]') _diag(final_G,exp(diagonal(final_G))) // CONVERTS logcholGstar TO cholGstar _lowertriangle(final_G) final_G = final_G*final_G' // CONVERTS cholGstar TO Gstar final_G = final_sigmaSquared:*final_G // CONVERTS Gstar to G // CONVERT PARAMETERS FROM G -> arctanG final_arctanG = final_G // OFF DIAGONAL TERMS for(row=1; row<=colsZ; row++) { for(col=row+1; col<=colsZ; col++) { final_arctanG[row,col] = atanh(final_arctanG[row,col]/sqrt(final_arctanG[row,row]*final_arctanG[col,col])) final_arctanG[col,row] = final_arctanG[row,col] } } // DIAGONAL TERMS _diag(final_arctanG,log(sqrt(diagonal(final_arctanG)))) // CONVERT R PARAMETERS TO (lnalpha,lntau,lnsigma) OR (lnphi,lnsigma) if (Rparameterization < 8) { // IOU MODEL // IOUstar -> (alpha, taustar) Rconversion(Rparameterization, optimizedParameters[1,numberOfGparas+1], optimizedParameters[1,numberOfGparas+2], /// 1, final_alpha=., final_tau=.) // CONVERTS taustar -> tau final_tau = final_tau*final_sigma final_Rparameters = (final_alpha,final_tau,final_sigmaSquared) // alpha, tau, sigma -> lnalpha, lntau, lnsigma optimizedParameters = (vech(final_arctanG)',ln(final_alpha),ln(final_tau),ln(final_sigma)) settings = (2,7,1) } else { // BROWNIAN-MOTION MODEL final_phi = optimizedParameters[1,numberOfGparas+1] final_phi = final_phi*final_sigmaSquared final_Rparameters = (final_phi,final_sigmaSquared) // phi, sigma -> lnphi, lnsigma optimizedParameters = (vech(final_arctanG)',ln(final_phi),ln(final_sigma)) settings = (2,9,1) } if (hasmissing(optimizedParameters)==1) { displayas("error") printf("One or more estimates of the variance parameters are invalid\n") exit(error(430)) } // CALCULATION OF COVARIANCE MATRIX FOR (vech(arctanG), ln R PARAMETERS) USING NR'S HESSIAN printf("\nComputing standard errors:\n") NR_REML_d2(0, optimizedParameters, X, Y, Z, panelInfo, timeVar, settings, sigmaSquared=., /// beta=., varbeta=., lnf=., gradient=., Hessian=.) final_beta = beta final_varbeta = varbeta // NEEDS TO BE CALCULATED BASED ON THE UNPROFILED REML // b AND V FOR (beta', vech(arctanG), lnalpha, lntau, lnsigma) b = (final_beta', optimizedParameters) cols_b = cols(b) numFE = rows(final_beta) V = J(cols_b, cols_b, 0) V[|1,1 \ numFE, numFE|] = final_varbeta // TABULATE THE RESULTS FOR THE FIXED EFFECTS tabulatingBeta(final_beta, final_varbeta, table_beta=.) // TEST IF THE Hessian MATRIX IS POSITIVE DEFINITE isPositiveDefinite = isPositiveDefinite(Hessian) if(isPositiveDefinite==-1 | isPositiveDefinite==0) { errorCode = 3353 table_theta = J(cols(optimizedParameters),3,.) table_theta = ((vech(final_G)\final_Rparameters'), table_theta) } else { V[|numFE+1, numFE+1 \ cols_b, cols_b|] = 2:*invsym(Hessian) // TABULATE THE RESULTS FOR THE VARIANCE PARAMETERS // CALCULATES THE UNTRANSFORMED PARAMETERS // USES THE DELTA METHOD TO CALCULATE THE STANDARD ERRORS OF THE UNTRANSFORMED PARAMETERS // CALCULATES THE CIs ON THE TRANSFORMED SCALE AND THEN BACK-TRANSFORMS tabTheta_sigmaSq(optimizedParameters, V[|numFE+1, numFE+1 \ cols_b, cols_b|], colsZ, numberOfGparas, /// table_theta=.) } st_rclear() st_numscalar("r(errorCode)", errorCode) st_numscalar("r(converged)", converged) st_numscalar("r(ll_reml)", ll_reml) st_numscalar("r(numberOfObs)", numberOfObs) st_numscalar("r(numberOfPanels)", numberOfPanels) st_numscalar("r(min_ni)", min_ni) st_numscalar("r(max_ni)", max_ni) st_matrix("r(b)", b) st_matrix("r(V)", V) st_matrix("r(table_beta)", table_beta) st_matrix("r(table_theta)", table_theta) st_matrix("r(G)", final_G) st_matrix("r(Rparameters)", final_Rparameters) } // END OF FUNCTION xtmixediou_estimation() void NR_REML_d2(real scalar todo, real rowvector parameters, /// real matrix X, real colvector Y, real matrix Z, /// real matrix panelInfo, real colvector timeVar, /// real rowvector settings, real scalar sigmaSquared, /// real colvector beta, real matrix CC_T, /// real scalar lnf, real rowvector gradient, real matrix Hessian) { real scalar panel, numberOfPanels, colsX, colsZ, numberOfParameters, deriv_r, deriv_s, ni, N real scalar numberOfGparas, numberOfWparas, numberOfRparas, minusLogDetR, lndet_augW0_part1, lndet_augWG_part1 real scalar row_r, col_r, row_s, col_s, l2_scalar, l3_scalar real scalar aLpha, tau, lnsigma, triangle, iter_algorithm real scalar Gparameterization, Rparameterization, maxMissing real rowvector Gparameters, Rparameters, g1, g2, g3, g2star, missing, Wparameters real colvector residuali, Yi, timeVari, l1 real matrix Xi, Zi, sw1_augW0_part4, Xstari, matrix_g1, matrix_g2 real matrix vech_sw1_augW0_part4, sum_sw1_augW0_part4, sw1_augWG_part4, WGzz, WGzr, WGxz real matrix cholG, invcholG, Gdot_r, Gdot_s, Gdotdot_rs real matrix invcholRi, invRi, Rdot_r, Rdot_s, Rdotdot_rs real matrix vech_H1, H1, H2, H3, M, invRZMZ_TinvR, MZ_TinvR, A1_r, A2_r, B, sum_H3rs real matrix vech_H2rs, Wi, timeVari_s, timeVari_t, oneMatrix, minMatrix, Gpositions, Rpositions, Vpositions pointer(real matrix) colvector p_matrix_H2r, p_vech_H3r, p_sum_H3r, p_Gdot, p_Gdotdot, p_invcholRi, p_Rdot, p_Rdotdot pointer(real matrix) colvector p_vech_H3rs pointer(function) rowvector fnarray_Gdot, fnarray_Gdotdot, fnarray_Rcov, fnarray_Rdot, fnarray_Rdotdot // EXTRACT FEATURES FROM THE DATA numberOfPanels = panelstats(panelInfo)[1] N = rows(X) colsX = cols(X) colsZ = cols(Z) numberOfGparas = 0.5*colsZ*(colsZ+1) // CALCULATED FROM THE TRIANGLE NUMBER SERIES // EXTRACT PARAMETERIZATION AND PARAMETERS Gparameterization = settings[1,1] Rparameterization = settings[1,2] iter_algorithm = settings[1,3] numberOfParameters = cols(parameters) numberOfRparas = numberOfParameters - numberOfGparas Gparameters = parameters[|1,1 \ 1,numberOfGparas|] Rparameters = parameters[|1,numberOfGparas+1 \ 1,numberOfParameters|] // EXTRACT cholG [CODE TESTED BY "xtIOU\test command\test cholGextraction for logcholGstar and arctanG parameterizations.do" cholG = invvech(Gparameters') if (Gparameterization==1) { // logcholGstar PARAMETERIZATION _diag(cholG,exp(diagonal(cholG))) // logcholGstar -> cholGstar; EXPONENTIATE THE DIAGONAL ELEMENTS _lowertriangle(cholG) lnsigma = 0 // i.e.; ln(1) = 0 } else { // arctanG PARAMETERIZATION // arctanG -> G // OFF-DIAGONAL ELEMENTS for(row_r=1; row_r<=colsZ; row_r++) { for(col_r=row_r+1; col_r<=colsZ; col_r++) { cholG[row_r,col_r] = tanh(cholG[row_r,col_r])*exp(cholG[row_r,row_r])*exp(cholG[col_r,col_r]) cholG[col_r,row_r] = cholG[row_r,col_r] } } // DIAGONAL ELEMENTS _diag(cholG,exp(2:*diagonal(cholG))) // G -> cholG cholG = cholesky(cholG) lnsigma = Rparameters[1,cols(Rparameters)] } invcholG = luinv(cholG) // MAY BE luinv(cholesky[(1/sigma^2):*G]) OR luinv(cholesky(G)) if (Rparameterization < 8) { // OBTAIN aLpha AND tau (OR taustar) Rconversion(Rparameterization, Rparameters[1,1], Rparameters[1,2], 1, aLpha=., tau=.) Wparameters = (aLpha, tau) } else if (Rparameterization == 8) { Wparameters = Rparameters[1,1] // phistar } else { Wparameters = exp(Rparameters[1,1]) // ln(phi) } numberOfWparas = cols(Wparameters) // INITIALIZE VECTORS AND MATRICES Gpositions = invvech(1::(0.5*numberOfGparas*(numberOfGparas+1))) Rpositions = invvech(1::(0.5*numberOfRparas*(numberOfRparas+1))) l1 = J(numberOfPanels,1,0) triangle = colsX+colsZ+1 triangle = 0.5*triangle*(triangle+1) vech_sw1_augW0_part4 = J(triangle,numberOfPanels,0) matrix_g1 = J(numberOfPanels,numberOfParameters,0) matrix_g2 = J(numberOfPanels,numberOfParameters,0) g3 = J(1,numberOfParameters,0) triangle = 0.5*numberOfParameters*(numberOfParameters+1) Vpositions = invvech(1::triangle) vech_H1 = J(triangle,numberOfPanels,0) vech_H2rs = J(triangle,numberOfPanels,0) H2 = J(numberOfParameters,numberOfParameters,0) H3 = J(numberOfParameters,numberOfParameters,0) // INITIALIZE POINTERS p_Gdot = J(numberOfGparas,1,NULL) p_Gdotdot = J(0.5*numberOfGparas*(numberOfGparas+1),1,NULL) p_invcholRi = J(numberOfPanels,1,NULL) p_sum_H2r = J(numberOfParameters, 1, NULL) p_sum_H3r = J(numberOfParameters, 1, NULL) // NON-NULL POINTERS (REQUIRES A SEPARATE ZERO MATRIX FOR EACH POINTER ELEMENT) triangle = 0.5*colsX*(colsX+1) p_matrix_H2r = &(J(colsX,numberOfPanels,0)) p_vech_H3r = &(J(triangle,numberOfPanels,0)) p_vech_H3rs = &(J(triangle,numberOfPanels,0)) for(deriv_r=2; deriv_r<=numberOfParameters; deriv_r++) { p_vech_H3r = p_vech_H3r \ &(J(triangle,numberOfPanels,0)) p_matrix_H2r = p_matrix_H2r \ &(J(colsX,numberOfPanels,0)) for(deriv_s=1; deriv_s<=deriv_r; deriv_s++) { p_vech_H3rs = p_vech_H3rs \ &(J(triangle,numberOfPanels,0)) } } // FUNCTION ARRAYS fnarray_Rcov = (&Browniancov(),&IOUcov()) fnarray_Gdot = (&Gdot_logcholG(), &Gdot_arctanG()) fnarray_Gdotdot = (&Gdotdot_logcholG(), &Gdotdot_arctanG()) fnarray_Rdot = (&Rdot_at(), &Rdot_ao(), &Rdot_et(), &Rdot_eo(), /// &Rdot_it(), &Rdot_io(), &Rdot_ek(), &Rdot_phi(), &Rdot_lnphi()) fnarray_Rdotdot = (&Rdotdot_at(), &Rdotdot_ao(), &Rdotdot_et(), &Rdotdot_eo(), /// &Rdotdot_it(), &Rdotdot_io(), &Rdotdot_ek(), &Rdotdot_phi(), &Rdotdot_lnphi()) // GENERATE Gdot AND Gdotdot col_r = 1 row_r = 0 for(deriv_r=1; deriv_r<=numberOfGparas; deriv_r++) { row_r = row_r + 1 if(row_r > colsZ) { col_r = col_r + 1 row_r = col_r } p_Gdot[deriv_r] = &((*fnarray_Gdot[Gparameterization])(Gparameters, row_r, col_r)) col_s = 1 row_s = 0 for(deriv_s=1; deriv_s<=deriv_r; deriv_s++) { row_s = row_s + 1 if(row_s > colsZ) { col_s = col_s + 1 row_s = col_s } p_Gdotdot[Gpositions[deriv_r,deriv_s]] = &((*fnarray_Gdotdot[Gparameterization])(Gparameters, row_r, col_r, row_s, col_s)) } // END OF deriv_s FOR-LOOP } // END OF deriv_r FOR-LOOP // CALCULATING lnf for(panel=1; panel<=numberOfPanels; panel++) { ni = panelInfo[panel,2] - panelInfo[panel,1] + 1 Zi = panelsubmatrix(Z,panel,panelInfo) Xi = panelsubmatrix(X,panel,panelInfo) Yi = panelsubmatrix(Y,panel,panelInfo) timeVari = panelsubmatrix(timeVar,panel,panelInfo) oneMatrix = J(ni,ni,1) minMatrix = makesymmetric(lowertriangle(timeVari' :* oneMatrix)) timeVari_s = timeVari :* oneMatrix timeVari_t = timeVari':* oneMatrix // GENERATE COVARIANCE MATRIX Wi (*fnarray_Rcov[numberOfWparas])(Wparameters, oneMatrix, timeVari_s, timeVari_t, minMatrix, Wi=.) p_invcholRi[panel] = &(luinv(cholesky(Wi + exp(2*lnsigma):*I(ni)))) // POINTER MUST POINT TO CALCULATION AND NOT MATRIX CALLED invcholRi //(i.e. IF p_invcholRi[panel]=invcholRi THEN A PROBLEM OCCURS: CONTENTS OF POINTER CHANGES AS invcholRi CHANGES invcholRi = *p_invcholRi[panel] invRi = cross(invcholRi,invcholRi) /********************************************************************* WOLFINGER LIKELIHOOD CALCULATIONS *********************************************************************/ // SECTION 3 OF "check - xtIOU" SHOWS THAT THE FOLLOWING GIVES -log|R| minusLogDetR = 2*quadcolsum(ln(diagonal(invcholRi))) // GENERATES MATRIX W AND THE LOG(DET(FIRST PARTITION OF AUGMENTED W0)) // STEP 1: GENERATE CROSS-PRODUCTS MATRIX W0 // STEP 2: GENERATE THE AUGMENTED MATRIX OF W0 // STEP 3: SWEEP THE AUGMENTED MATRIX // W = C-B'A^-1B OF THE SWEEP; THE FOURTH PARTITION OF THE SWEPT MATRIX crossProductMatrices(Xi, colsX, Yi, Zi, colsZ, cholG, invcholRi, lndet_augW0_part1=., sw1_augW0_part4=.) // l1 = ln|I+L'*Z'*invR*Z*L| + log|R| // = log|V| - log|R| + log|R| l1[panel] = lndet_augW0_part1 - minusLogDetR // EXITS EARLY IF ONE PANEL HAS A MISSING VALUE maxMissing = hasmissing(l1[panel]) if (maxMissing==1) break // TO CALCULATE l2 AND l3 YOU MUST SUM sw1_augW0_part4 OVER THE m SUBJECTS // e.g. SEE THE FORMULAE GIVEN TO Lindstrom and Bates 1988 // THE MATRIX IS SYMMETRIC vech_sw1_augW0_part4[,panel] = vech(sw1_augW0_part4) } // END OF panel FOR-LOOP sum_sw1_augW0_part4 = invvech(quadrowsum(vech_sw1_augW0_part4)) // SWEEP THE APPROPRIATE SUBMATRIX OF sw1_augW0_part4 TO OBTAIN l2 AND l3 // APPLIED TO THE sum OF SUBMATRIX OF W0 sweep_partitionK(sum_sw1_augW0_part4[|1,1\colsX,colsX|], sum_sw1_augW0_part4[|1,1+colsX+colsZ \colsX,1+colsX+colsZ|], sum_sw1_augW0_part4[1+colsX+colsZ,1+colsX+colsZ], /// CC_T=., beta=., l2_scalar=., l3_scalar=.) sigmaSquared = l2_scalar/(N-colsX) if(maxMissing==0 & l2_scalar !=. & l3_scalar !=.) { /********************************************************************* WOLFINGER DERIVATIVE CALCULATIONS *********************************************************************/ for(panel=1; panel<=numberOfPanels; panel++) { Zi = panelsubmatrix(Z,panel,panelInfo) Xi = panelsubmatrix(X,panel,panelInfo) Yi = panelsubmatrix(Y,panel,panelInfo) timeVari = panelsubmatrix(timeVar,panel,panelInfo) invcholRi = *p_invcholRi[panel] invRi = cross(invcholRi,invcholRi) ni = rows(Yi) oneMatrix = J(ni,ni,1) minMatrix = makesymmetric(lowertriangle(timeVari' :* oneMatrix)) timeVari_s = timeVari :* oneMatrix timeVari_t = timeVari':* oneMatrix // DERIVATIVES OF R (*fnarray_Rdot[Rparameterization])(ni, Rparameters, timeVari_s, timeVari_t, minMatrix, p_Rdot=NULL) (*fnarray_Rdotdot[Rparameterization])(ni, Rparameters, timeVari_s, timeVari_t, minMatrix, p_Rdotdot=NULL) Xstari = Xi*cholesky(CC_T) // SAME AS POPULATION LEVELS residuali = Yi-Xi*beta // DERIVATIVE CALCULATIONS USING WOLFINGER 1994 // GENERATE MATRIX WG - this matrix is required for the second derivatives too. crossProductMatrices(Xstari, colsX, residuali, Zi, colsZ, cholG, invcholRi, lndet_augWG_part1=., sw1_augWG_part4=.) // SUBMATRICES OF WG WGzz = sw1_augWG_part4[|colsX+1,colsX+1 \ colsX+colsZ,colsX+colsZ|] WGzr = sw1_augWG_part4[|colsX+1,colsX+colsZ+1 \ colsX+colsZ,colsX+colsZ+1|] WGxz = sw1_augWG_part4[|1,colsX+1 \ colsX,colsX+colsZ|] // DERIVATIVES W.R.T. PARAMETERS OF MATRIX G for(deriv_r=1; deriv_r<=numberOfGparas; deriv_r++) { Gdot_r = *p_Gdot[deriv_r] matrix_g1[panel,deriv_r] = trace(WGzz,Gdot_r) matrix_g2[panel,deriv_r] = - WGzr'*Gdot_r*WGzr (*p_matrix_H2r[deriv_r])[,panel] = WGxz*Gdot_r*WGzr (*p_vech_H3r[deriv_r])[,panel] = vech(WGxz*Gdot_r*WGxz') for(deriv_s=1; deriv_s<=deriv_r; deriv_s++) { Gdot_s = *p_Gdot[deriv_s] Gdotdot_rs = *p_Gdotdot[Gpositions[deriv_r,deriv_s]] vech_H2rs[Vpositions[deriv_r,deriv_s],panel] = 2*WGzr'*Gdot_r*WGzz*Gdot_s*WGzr - WGzr'*Gdotdot_rs*WGzr (*p_vech_H3rs[Vpositions[deriv_r,deriv_s]])[,panel] = vech(2*WGxz*Gdot_r*WGzz*Gdot_s*WGxz' - WGxz*Gdotdot_rs*WGxz') // MUST FILL IN THE LOWER TRIANGLE AS makesymmetric REFLECTS THE BELOW DIAGONAL TO ABOVE // SECOND DERIVATIVES OF AN UNSTRUCTURED G MATRIX ALWAYS RESULTS IN A ZERO VECTOR vech_H1[Vpositions[deriv_r,deriv_s],panel] = -trace(WGzz*Gdot_r,WGzz*Gdot_s) + trace(WGzz,Gdotdot_rs) } // END OF deriv_s FOR-LOOP } // END OF deriv_r FOR-LOOP // DERIVATIVES W.R.T. PARAMETERS OF COVARIANCE MATRIX R M = luinv(cross(invcholG,invcholG) + Zi'*invRi*Zi) invRZMZ_TinvR = invRi*Zi*M*Zi'*invRi MZ_TinvR = M*Zi'*invRi for(deriv_r=numberOfGparas+1; deriv_r<=numberOfParameters; deriv_r++) { col_r = deriv_r-numberOfGparas Rdot_r = *p_Rdot[col_r] (*p_matrix_H2r[deriv_r])[,panel] = Xstari'*invRi*Rdot_r*invRi*residuali /// -Xstari'*invRZMZ_TinvR*Rdot_r*invRi*residuali /// -Xstari'*invRi*Rdot_r*invRZMZ_TinvR*residuali /// +Xstari'*invRZMZ_TinvR*Rdot_r*invRZMZ_TinvR*residuali (*p_vech_H3r[deriv_r])[,panel] = vech(Xstari'*invRi*Rdot_r*invRi*Xstari /// -Xstari'*invRZMZ_TinvR*Rdot_r*invRi*Xstari /// -Xstari'*invRi*Rdot_r*invRZMZ_TinvR*Xstari /// +Xstari'*invRZMZ_TinvR*Rdot_r*invRZMZ_TinvR*Xstari) // GRADIENT COMPONENTS CORRESPONDING TO PARAMETERS OF R matrix_g1[panel,deriv_r] = trace(invRi*Rdot_r) - /// trace(MZ_TinvR*Rdot_r*invRi*Zi) matrix_g2[panel,deriv_r] = -residuali'*invRi*Rdot_r*invRi*residuali /// +2:*residuali'*invRZMZ_TinvR*Rdot_r*invRi*residuali /// -residuali'*invRZMZ_TinvR*Rdot_r*invRZMZ_TinvR*residuali for(deriv_s=numberOfGparas+1; deriv_s<=deriv_r; deriv_s++) { col_s = deriv_s-numberOfGparas Rdot_s = *p_Rdot[col_s] Rdotdot_rs = *p_Rdotdot[Rpositions[col_r,col_s]] // SIMPLIFYING COMPONENTS vech_H2rs[Vpositions[deriv_r,deriv_s],panel] = 2:*residuali'*invRi*Rdot_r*invRi*Rdot_s*invRi*residuali /// -2:*residuali'*invRZMZ_TinvR*Rdot_r*invRi*Rdot_s*invRi*residuali /// -2:*residuali'*invRi*Rdot_r*invRZMZ_TinvR*Rdot_s*invRi*residuali /// +2:*residuali'*invRZMZ_TinvR*Rdot_r*invRZMZ_TinvR*Rdot_s*invRi*residuali /// -2:*residuali'*invRi*Rdot_r*invRi*Rdot_s*invRZMZ_TinvR*residuali /// +2:*residuali'*invRZMZ_TinvR*Rdot_r*invRi*Rdot_s*invRZMZ_TinvR*residuali /// +2:*residuali'*invRi*Rdot_r*invRZMZ_TinvR*Rdot_s*invRZMZ_TinvR*residuali /// -2:*residuali'*invRZMZ_TinvR*Rdot_r*invRZMZ_TinvR*Rdot_s*invRZMZ_TinvR*residuali /// -residuali'*invRi*Rdotdot_rs*invRi*residuali /// +2:*residuali'*invRi*Rdotdot_rs*invRZMZ_TinvR*residuali /// -residuali'*invRZMZ_TinvR*Rdotdot_rs*invRZMZ_TinvR*residuali (*p_vech_H3rs[Vpositions[deriv_r,deriv_s]])[,panel] = vech(2:*Xstari'*invRi*Rdot_r*invRi*Rdot_s*invRi*Xstari /// -2:*Xstari'*invRZMZ_TinvR*Rdot_r*invRi*Rdot_s*invRi*Xstari /// -2:*Xstari'*invRi*Rdot_r*invRZMZ_TinvR*Rdot_s*invRi*Xstari /// +2:*Xstari'*invRZMZ_TinvR*Rdot_r*invRZMZ_TinvR*Rdot_s*invRi*Xstari /// -2:*Xstari'*invRi*Rdot_r*invRi*Rdot_s*invRZMZ_TinvR*Xstari /// +2:*Xstari'*invRZMZ_TinvR*Rdot_r*invRi*Rdot_s*invRZMZ_TinvR*Xstari /// +2:*Xstari'*invRi*Rdot_r*invRZMZ_TinvR*Rdot_s*invRZMZ_TinvR*Xstari /// -2:*Xstari'*invRZMZ_TinvR*Rdot_r*invRZMZ_TinvR*Rdot_s*invRZMZ_TinvR*Xstari /// -Xstari'*invRi*Rdotdot_rs*invRi*Xstari /// +2:*Xstari'*invRi*Rdotdot_rs*invRZMZ_TinvR*Xstari /// -Xstari'*invRZMZ_TinvR*Rdotdot_rs*invRZMZ_TinvR*Xstari) // FILLING IN MATRIX H vech_H1[Vpositions[deriv_r,deriv_s],panel] = - trace(invRi*Rdot_r*invRi*Rdot_s) /// + trace(MZ_TinvR*Rdot_r*invRi*Rdot_s*invRi*Zi) /// + trace(MZ_TinvR*Rdot_s*invRi*Rdot_r*invRi*Zi) /// - trace(MZ_TinvR*Rdot_r*invRZMZ_TinvR*Rdot_s*invRi*Zi) /// + trace(invRi*Rdotdot_rs) /// - trace(MZ_TinvR*Rdotdot_rs*invRi*Zi) } // END OF deriv_s FOR-LOOP } // END OF deriv_r FOR-LOOP // CROSS DERIVATIVES d/dG(d/dR) // e.g. H1[r,s] DERIVATIVE OF L1 W.R.T. PARAMETER r FROM R AND PARAMETER s FROM G for(deriv_r=numberOfGparas+1; deriv_r<=numberOfGparas+numberOfRparas; deriv_r++) { Rdot_r = *p_Rdot[deriv_r-numberOfGparas] for(deriv_s=1; deriv_s<=numberOfGparas; deriv_s++) { Gdot_s = *p_Gdot[deriv_s,1] // SIMPLYFYING COMPONENTS A1_r = residuali'*invRi*Rdot_r*invRi*Zi /// -residuali'*invRZMZ_TinvR*Rdot_r*invRi*Zi A2_r = Xstari'*invRi*Rdot_r*invRi*Zi /// -Xstari'*invRZMZ_TinvR*Rdot_r*invRi*Zi B = MZ_TinvR*Zi - I(colsZ) vech_H2rs[Vpositions[deriv_r,deriv_s],panel] = 2*A1_r*B*Gdot_s*B'*Zi'*invRi*residuali (*p_vech_H3rs[Vpositions[deriv_r,deriv_s]])[,panel] = vech(2*A2_r*B*Gdot_s*B'*Zi'*invRi*Xstari) // CROSS DERIVATIVES vech_H1[Vpositions[deriv_r,deriv_s],panel] = - trace(Zi'*invRi*Rdot_r*invRi*Zi*Gdot_s) /// + 2*trace(Zi'*invRi*Rdot_r*invRi*Zi*Gdot_s*Zi'*invRi*Zi*M) /// - trace(Zi'*invRi*Rdot_r*invRZMZ_TinvR*Zi*Gdot_s*Zi'*invRi*Zi*M) } // END OF deriv_s FOR-LOOP } // END OF deriv_r FOR-LOOP // EXITS EARLY IF ONE PANEL HAS A MISSING VALUE missing = hasmissing(matrix_g1[panel,]) missing = missing, hasmissing(matrix_g2[panel,]) missing = missing, hasmissing(vech_H1[,panel]) missing = missing, hasmissing(vech_H2rs[,panel]) maxMissing = max(missing) if (maxMissing==1) break } // END OF panel FOR-LOOP if(maxMissing==0) { H1 = invvech(quadrowsum(vech_H1)) vech_H2rs = quadrowsum(vech_H2rs) // AVOIDS IF STATEMENTS IN THE PREVIOUS LOOP AND COLLECTS RESULTS BY LOOPING THROUGH THE PARAMETERS AGAIN for(deriv_r=1; deriv_r<=numberOfParameters; deriv_r++) { p_sum_H2r[deriv_r] = &(quadrowsum(*p_matrix_H2r[deriv_r])) p_sum_H3r[deriv_r] = &(invvech(quadrowsum(*p_vech_H3r[deriv_r]))) g3[deriv_r] = -trace(*p_sum_H3r[deriv_r]) for(deriv_s=1; deriv_s<=deriv_r; deriv_s++) { sum_H3rs = invvech(quadrowsum(*p_vech_H3rs[Vpositions[deriv_r,deriv_s]])) H2[deriv_r,deriv_s] = vech_H2rs[Vpositions[deriv_r,deriv_s]] - 2*(*p_sum_H2r[deriv_r]')*(*p_sum_H2r[deriv_s]) H3[deriv_r,deriv_s] = trace(sum_H3rs - (*p_sum_H3r[deriv_r])*(*p_sum_H3r[deriv_s])) } // END OF deriv_s FOR-LOOP } // END OF deriv_r FOR-LOOP _makesymmetric(H1) // SAVES FILLING IN THE OTHER ENTRIES _makesymmetric(H2) // SAVES FILLING IN THE OTHER ENTRIES _makesymmetric(H3) // SAVES FILLING IN THE OTHER ENTRIES g1 = quadcolsum(matrix_g1) g2 = quadcolsum(matrix_g2) // CALCULATES lnf, gradient AND THE Hessian MATRIX if(Gparameterization == 1) { // logcholGstar AND taustar OR omegastar lnf = quadcolsum(l1) + (N-colsX)*log(l2_scalar) + l3_scalar /// + (N-colsX) + (N-colsX)*ln(2*pi()/(N-colsX)) g2star = (1/sigmaSquared):*g2 gradient = g1 + g2star + g3 if(iter_algorithm==1) { // NR Hessian = H1 + (1/sigmaSquared):*H2 - (1/(N-colsX)):*g2star'*g2star + H3 } else if(iter_algorithm==2) { // FS Hessian = -H1 + H3 - (1/(N-colsX)):*g2star'*g2star } else { // AI Hessian = 0.5*(1/sigmaSquared):*H2 + H3 - (1/(N-colsX)):*g2star'*g2star } } else { // arctanG AND tau OR omega AND sigma lnf = quadcolsum(l1) + l2_scalar + l3_scalar + (N-colsX)*ln(2*pi()) gradient = g1 + g2 + g3 if(iter_algorithm==1) { // NR Hessian = H1 + H2 + H3 } else if(iter_algorithm==2) { // FS Hessian = -H1 + H3 } else { // AI Hessian = 0.5:*H2 + H3 } } } // END OF if g1,g2,vech_H1,vech_H2rs HAS A MISSING VALUE else { lnf = . gradient = . Hessian = . } } // END OF if(maxMissing==0 & l2_scalar !=. & l3_scalar !=.) else { lnf = . gradient = . Hessian = . } } // END OF FUNCTION NR_REML_d2() void IOUcov(real rowvector Wparameters, real matrix oneMatrix, real matrix timeVari_s, real matrix timeVari_t, real matrix minMatrix, real matrix Wi) { real scalar aLpha, tau aLpha = Wparameters[1,1] tau = Wparameters[1,2] Wi = (tau^2)/((aLpha^3)*2):*(2*aLpha:*minMatrix + exp(-aLpha:*timeVari_s) + exp(-aLpha:*timeVari_t) - oneMatrix /// - exp(-aLpha:*abs(timeVari_t-timeVari_s))) } // END OF FUNCTION IOUcov() void Browniancov(real rowvector Wparameters, real matrix oneMatrix, real matrix timeVari_s, real matrix timeVari_t, real matrix minMatrix, real matrix Wi) { real scalar phi phi = Wparameters[1,1] Wi = phi:*minMatrix } // END OF FUNCTION Browniancov() // FUNCTION PERFORMS THE FOLLOWING: // GENERATES CROSS-PRODUCT MATRIX W0 // GENERATES THE AUGMENTATION OF W0 // SWEEPS THE AUGMENTATION OF W0 ON THE FIRST PARTITION void crossProductMatrices(real matrix A, real scalar colsA, real colvector B, real matrix C, real scalar colsC, real matrix cholD, /// real matrix invcholE, real scalar lndet_augW_part1, real matrix sw1_augW_part4) { real matrix invcholEA, invcholEC, invcholEB, W, augW_part1, augW_part2, sw1_augW_part1, sw1_augW_part2 invcholEA = invcholE*A // invcholE*A invcholEC = invcholE*C // invcholE*C invcholEB = invcholE*B // invcholE*B // GENERATE CROSS-PRODUCT MATRIX W W = J(1+colsA+colsC,1+colsA+colsC,0) W[|1,1\colsA,colsA|] = cross(invcholEA,invcholEA) // A'*invcholE'*invcholE*A W[|1+colsA,1\colsA+colsC,colsA|] = cross(invcholEC,invcholEA) // C'*invcholE'*invcholE*A W[|1+colsA,1+colsA\colsA+colsC,colsA+colsC|] = cross(invcholEC,invcholEC) // C'invcholE'*invcholE*C W[|1+colsA+colsC,1\1+colsA+colsC,colsA|] = cross(invcholEB,invcholEA) // B'*invcholE'*invcholE*A W[|1+colsA+colsC,1+colsA\1+colsA+colsC,colsA+colsC|] = cross(invcholEB,invcholEC) // B'*invcholE'*invcholE*C W[1+colsA+colsC,1+colsA+colsC] = cross(invcholEB,invcholEB) // B'*invcholE'*invcholE*B _makesymmetric(W) // SAVES FILLING IN THE OTHER ENTRIES // GENERATE AUGMENTATION OF W - IN 4 PARTITIONS // CODING CHECKED - SEE SECTION 2 OF "check - xtIOU.do" // aug W = [augW_part1 augW_part2 // augW_part3 augW_part4] // ONLY REQUIRE THE UPPER TRIANGLE (PARTS 1, 3 AND 4) // augW_part4 = W augW_part1 = I(colsC) + cross(invcholEC*cholD,invcholEC*cholD) augW_part2 = cholD'*W[|colsA+1,1 \ colsA+colsC,colsA+colsC+1|] // SWEEP THE AUGMENTED W MATRIX AND REPORT THE SUM THE LOG OF THE POSITIVE PIVOTS DURING THE SWEEP sweep_partitionK(augW_part1, augW_part2, W, sw1_augW_part1=., sw1_augW_part2=., sw1_augW_part4=., lndet_augW_part1=.) } // END OF FUNCTION crossProductMatrices void sweep_partitionK(real matrix partitionKK, real matrix partitionKJ, real matrix partitionJJ, real matrix sweptKK, real matrix sweptKJ, /// real matrix sweptJJ, real scalar lndet_partitionKK) { sweptKK = luinv(partitionKK) sweptKJ = sweptKK*partitionKJ sweptJJ = partitionJJ - cross(partitionKJ,sweptKJ) // partitionJJ - partitionJK*sweptKK*partitionKJ // SUM THE LOG OF THE POSITIVE PIVOTS DURING THE SWEEP lndet_partitionKK = ln(det(partitionKK)) } // END OF FUNCTION sweep_partitionK real matrix Gdot_logcholG(real rowvector Gparameters, real scalar row, real scalar col) { real scalar colsZ real matrix logcholG, Gdot logcholG = invvech(Gparameters') colsZ = cols(logcholG) Gdot = J(colsZ,colsZ,0) if(row==col) { Gdot[|row,row \ row,colsZ|] = exp(logcholG[row,row]):*logcholG[|row,row \ colsZ,row|]' Gdot[row,row] = exp(2*logcholG[row,row]) } else { Gdot[|row,col \ row,colsZ|] = logcholG[|col,col \ colsZ,col|]' Gdot[row,col] = exp(logcholG[col,col]) } Gdot[row, row] = 2*Gdot[row, row] Gdot[|1,row \ colsZ,row|] = Gdot[|row,1 \ row,colsZ|]' return(Gdot) } // END OF FUNCTION Gdot_logcholG real matrix Gdotdot_logcholG(real rowvector Gparameters, real scalar row_r, real scalar col_r, real scalar row_s, real scalar col_s) { real scalar colsZ real matrix logcholG, Gdotdot logcholG = invvech(Gparameters') colsZ = cols(logcholG) Gdotdot = J(colsZ,colsZ,0) // deriv_r==deriv_s if(row_r==row_s & col_r==col_s) { if(row_r==col_r) { Gdotdot[|row_r,row_r \ colsZ,row_r|] = exp(logcholG[row_r,row_r]):*logcholG[|row_r,row_r \ colsZ,row_r|] Gdotdot[|row_r,row_r \ row_r,colsZ|] = Gdotdot[|row_r,row_r \colsZ,row_r|]' Gdotdot[row_r,row_r] = 4*exp(2*logcholG[row_r,row_r]) } else { Gdotdot[row_r,row_r] = 2 } } else if(row_r!=row_s & col_r==col_s) { // e.g. row_r=3,col_r=2 and row_s=2,col_s=2 if(row_s==col_s) { Gdotdot[row_r,col_r] = exp(logcholG[row_s,row_s]) Gdotdot[col_r,row_r] = Gdotdot[row_r,col_r] } else { Gdotdot[row_r,row_s] = 1 Gdotdot[row_s,row_r] = 1 } } // ELSE A ZERO MATRIX (i.e. row_r!=row_s & col_r!=col_s) return(Gdotdot) } // END OF FUNCTION Gdotdot_logcholG real matrix Gdot_arctanG(real rowvector Gparameters, real scalar row, real scalar col) { real scalar colsZ, index real matrix arctanG, Gdot arctanG = invvech(Gparameters') colsZ = cols(arctanG) Gdot = J(colsZ,colsZ,0) if(row==col) { for(index=1; index<=colsZ; index++) { Gdot[row,index] = exp(arctanG[index,index]+arctanG[row,row])*tanh(arctanG[row,index]) } Gdot[row,row] = 2*exp(2*arctanG[row,row]) Gdot[|1,row \ colsZ,row|] = Gdot[|row,1 \ row,colsZ|]' } else { Gdot[row,col] = exp(arctanG[row,row]+arctanG[col,col])*(cosh(arctanG[row,col])^-2) // sech[x] = 1/cosh[x] Gdot[col,row] = Gdot[row,col] } return(Gdot) } // END OF FUNCTION Gdot_arctanG real matrix Gdotdot_arctanG(real rowvector Gparameters, real scalar row_r, real scalar col_r, real scalar row_s, real scalar col_s) { real scalar colsZ, index real matrix arctanG, Gdotdot arctanG = invvech(Gparameters') colsZ = cols(arctanG) Gdotdot = J(colsZ,colsZ,0) if(row_r==row_s & col_r==col_s) { // deriv_r==deriv_s if(row_r==col_r) { // OPTION 1 for(index=1; index<=colsZ; index++) { Gdotdot[row_r,index] = exp(arctanG[index,index]+arctanG[row_r,row_r])*tanh(arctanG[row_r,index]) } Gdotdot[row_r,row_r] = 4*exp(2*arctanG[row_r,row_r]) Gdotdot[|1,row_r \ colsZ,row_r|] = Gdotdot[|row_r,1 \ row_r,colsZ|]' } else { // OPTION 2 Gdotdot[row_r,col_r] = -2*exp(arctanG[row_r,row_r]+arctanG[col_r,col_r])*(1/cosh(arctanG[row_r,col_r]))^2*tanh(arctanG[row_r,col_r]) Gdotdot[col_r,row_r] = Gdotdot[row_r,col_r] } } else { // deriv_r!=deriv_s if(col_r==row_s & row_s==col_s) { Gdotdot[row_r,col_r] = exp(arctanG[row_r,row_r]+arctanG[col_r,col_r])*(1/cosh(arctanG[row_r,col_r]))^2 Gdotdot[col_r,row_r] = Gdotdot[row_r,col_r] } else if(col_r==row_s & row_r==col_r) { Gdotdot[row_r,col_s] = exp(arctanG[row_r,row_r]+arctanG[col_s,col_s])*(1/cosh(arctanG[row_r,col_s]))^2 Gdotdot[col_s,row_r] = Gdotdot[row_r,col_s] } else if(row_r==col_r & row_s==col_s) { // OUTER else IMPLIES THAT row_r!=row_s AND col_r!=col_s Gdotdot[row_r,row_s] = exp(arctanG[row_r,row_r]+arctanG[row_s,row_s])*tanh(arctanG[row_r,row_s]) Gdotdot[row_s,row_r] = Gdotdot[row_r,row_s] } } // ELSE ZERO MATRIX return(Gdotdot) } // END OF FUNCTION Gdotdot_arctanG void Rdot_at(real scalar ni, real rowvector Rparameters, /// real matrix timeVari_s, real matrix timeVari_t, real matrix minMatrix, /// pointer(real matrix) colvector p_Rdot) { real scalar aLpha, tau real matrix oneMatrix aLpha = Rparameters[1] tau = Rparameters[2] oneMatrix = J(ni,ni,1) p_Rdot = J(2,1,NULL) // DIFFERENTIATE WITH RESPECT TO aLpha p_Rdot[1] = &(makesymmetric( /// (tau^2)/(2*(aLpha^4)):*(3:*oneMatrix - 3:*exp(-aLpha:*timeVari_s) - 3:*exp(-aLpha:*timeVari_t) + 3:*exp(-aLpha:*abs(timeVari_t-timeVari_s)) /// - exp(-aLpha:*timeVari_s):*(aLpha:*timeVari_s) - exp(-aLpha:*timeVari_t):*(aLpha:*timeVari_t) /// + exp(-aLpha:*abs(timeVari_t-timeVari_s)):*(aLpha:*abs(timeVari_t-timeVari_s)) - 4*aLpha:*minMatrix) /// )) // DIFFERENTIATE WITH RESPECT TO tau p_Rdot[2] = &(makesymmetric( /// tau/(aLpha^3):*(2*aLpha:*minMatrix + exp(-aLpha:*timeVari_s) + exp(-aLpha:*timeVari_t) - /// oneMatrix - exp(-aLpha:*abs(timeVari_t-timeVari_s))) /// )) } // END OF FUNCTION Rdot_at() void Rdotdot_at(real scalar ni, real rowvector Rparameters, /// real matrix timeVari_s, real matrix timeVari_t, real matrix minMatrix, /// pointer(real matrix) colvector p_Rdotdot) { real scalar aLpha, tau real matrix oneMatrix aLpha = Rparameters[1] tau = Rparameters[2] oneMatrix = J(ni,ni,1) p_Rdotdot = J(3,1,NULL) // DIFFERENTIATE WITH RESPECT TO aLpha AND aLpha p_Rdotdot[1] = &(makesymmetric( /// (tau^2)/(2*(aLpha^5)):*( /// (aLpha^2):*(exp(-aLpha:*timeVari_s):*timeVari_s:*timeVari_s + exp(-aLpha:*timeVari_t):*timeVari_t:*timeVari_t /// - exp(-aLpha:*abs(timeVari_t-timeVari_s)):*abs(timeVari_t-timeVari_s):*abs(timeVari_t-timeVari_s)) /// - 6*aLpha:*(-exp(-aLpha:*timeVari_s):*timeVari_s - exp(-aLpha:*timeVari_t):*timeVari_t /// + exp(-aLpha:*abs(timeVari_t-timeVari_s)):*abs(timeVari_t-timeVari_s) + 2:*minMatrix) /// + 12:*(-1:*oneMatrix + exp(-aLpha:*timeVari_s) + exp(-aLpha:*timeVari_t) - exp(-aLpha:*abs(timeVari_t-timeVari_s)) + 2*aLpha:*minMatrix)) /// )) // DIFFERENTIATE WITH RESPECT TO aLpha AND tau p_Rdotdot[2] = &(makesymmetric( /// tau/(aLpha^4):*(3:*oneMatrix - 3:*exp(-aLpha:*timeVari_s) - 3:*exp(-aLpha:*timeVari_t) + 3:*exp(-aLpha:*abs(timeVari_t-timeVari_s)) /// - exp(-aLpha:*timeVari_s):*aLpha:*timeVari_s - exp(-aLpha:*timeVari_t):*aLpha:*timeVari_t /// + exp(-aLpha:*abs(timeVari_t-timeVari_s)):*aLpha:*abs(timeVari_t-timeVari_s) - 4*aLpha:*minMatrix) /// )) // DIFFERENTIATE WITH RESPECT TO tau AND tau p_Rdotdot[3] = &(makesymmetric( /// 1/(aLpha^3):*(2*aLpha:*minMatrix + exp(-aLpha:*timeVari_s) + exp(-aLpha:*timeVari_t) /// - oneMatrix - exp(-aLpha:*abs(timeVari_t-timeVari_s))) /// )) } // END OF FUNCTION Rdotdot_at() void Rdot_ao(real scalar ni, real rowvector Rparameters, /// real matrix timeVari_s, real matrix timeVari_t, real matrix minMatrix, /// pointer(real matrix) colvector p_Rdot) { real scalar aLpha, omega real matrix oneMatrix aLpha = Rparameters[1] omega = Rparameters[2] oneMatrix = J(ni,ni,1) p_Rdot = J(2,1,NULL) // DIFFERENTIATE WITH RESPECT TO aLpha p_Rdot[1] = &(makesymmetric( /// (omega/(2*aLpha^2)):*(oneMatrix + exp(-aLpha:*abs(timeVari_t-timeVari_s)) - exp(-aLpha:*timeVari_s):*(oneMatrix + aLpha:*timeVari_s) /// - exp(-aLpha:*timeVari_t):*(oneMatrix + aLpha:*timeVari_t) + exp(-aLpha:*abs(timeVari_t-timeVari_s)):*aLpha:*abs(timeVari_t-timeVari_s)) /// )) // DIFFERENTIATE WITH RESPECT TO omega p_Rdot[2] = &(makesymmetric( /// (1/(2*aLpha)):*(2*aLpha:*minMatrix + exp(-aLpha:*timeVari_s) + exp(-aLpha:*timeVari_t) - oneMatrix - exp(-aLpha:*abs(timeVari_t-timeVari_s))) )) } // END OF FUNCTION Rdot_ao() void Rdotdot_ao(real scalar ni, real rowvector Rparameters, /// real matrix timeVari_s, real matrix timeVari_t, real matrix minMatrix, /// pointer(real matrix) colvector p_Rdotdot) { real scalar aLpha, omega real matrix oneMatrix aLpha = Rparameters[1] omega = Rparameters[2] oneMatrix = J(ni,ni,1) p_Rdotdot = J(3,1,NULL) // DIFFERENTIATE WITH RESPECT TO aLpha AND aLpha p_Rdotdot[1] = &(makesymmetric( /// (-omega/(2*aLpha^3)):*exp(-aLpha:*(timeVari_s+timeVari_t+abs(timeVari_t-timeVari_s))):* /// (2*exp(aLpha:*(timeVari_s+timeVari_t)) + 2*exp(aLpha:*(timeVari_s+timeVari_t+abs(timeVari_t-timeVari_s))) /// - exp(aLpha:*(timeVari_t+abs(timeVari_t-timeVari_s))):*(2:*oneMatrix+aLpha:*timeVari_s:*(2:*oneMatrix+aLpha:*timeVari_s)) /// - exp(aLpha:*(timeVari_s+abs(timeVari_t-timeVari_s))):*(2:*oneMatrix+aLpha:*timeVari_t:*(2:*oneMatrix+aLpha:*timeVari_t)) /// + 2*exp(aLpha:*(timeVari_s+timeVari_t)):*aLpha:*abs(timeVari_t-timeVari_s) /// + exp(aLpha:*(timeVari_s+timeVari_t)):*(aLpha^2):*(abs(timeVari_t-timeVari_s):*abs(timeVari_t-timeVari_s))) /// )) // DIFFERENTIATE WITH RESPECT TO aLpha AND omega p_Rdotdot[2] = &(makesymmetric( /// (1/(2*aLpha^2)):*(oneMatrix + exp(-aLpha:*abs(timeVari_t-timeVari_s)) - exp(-aLpha:*timeVari_s):*(oneMatrix+aLpha:*timeVari_s) /// - exp(-aLpha:*timeVari_t):*(oneMatrix+aLpha:*timeVari_t) + exp(-aLpha:*abs(timeVari_t-timeVari_s)):*aLpha:*abs(timeVari_t-timeVari_s)) /// )) // DIFFERENTIATE WITH RESPECT TO omega AND omega p_Rdotdot[3] = &( /// J(ni,ni,0) ) } // END OF FUNCTION Rdotdot_ao() void Rdot_et(real scalar ni, real rowvector Rparameters, /// real matrix timeVari_s, real matrix timeVari_t, real matrix minMatrix, /// pointer(real matrix) colvector p_Rdot) { real scalar aLpha, tau real matrix oneMatrix aLpha = exp(Rparameters[1]) tau = Rparameters[2] oneMatrix = J(ni,ni,1) p_Rdot = J(2,1,NULL) // DIFFERENTIATE WITH RESPECT TO eta p_Rdot[1] = &(makesymmetric( /// ((tau^2)/(2*(aLpha^3))):*(3:*(oneMatrix - exp(-aLpha:*timeVari_s) - exp(-aLpha:*timeVari_t) + exp(-aLpha:*abs(timeVari_t-timeVari_s))) /// + aLpha:*(-exp(-aLpha:*timeVari_s):*timeVari_s - exp(-aLpha:*timeVari_t):*timeVari_t + exp(-aLpha:*abs(timeVari_t-timeVari_s)):*abs(timeVari_t-timeVari_s) - 4:*minMatrix)) )) // DIFFERENTIATE WITH RESPECT TO tau p_Rdot[2] = &(makesymmetric( /// (tau/aLpha^3):*(2*aLpha:*minMatrix + exp(-aLpha:*timeVari_s) + exp(-aLpha:*timeVari_t) - oneMatrix - exp(-aLpha:*abs(timeVari_t-timeVari_s))) )) } // END OF FUNCTION Rdot_et() void Rdotdot_et(real scalar ni, real rowvector Rparameters, /// real matrix timeVari_s, real matrix timeVari_t, real matrix minMatrix, /// pointer(real matrix) colvector p_Rdotdot) { real scalar aLpha, tau real matrix oneMatrix aLpha = exp(Rparameters[1]) tau = Rparameters[2] oneMatrix = J(ni,ni,1) p_Rdotdot = J(3,1,NULL) // DIFFERENTIATE WITH RESPECT TO eta AND eta p_Rdotdot[1] = &(makesymmetric( /// ((tau^2)/(2*(aLpha)^3)):*(9:*(-oneMatrix + exp(-aLpha:*timeVari_s) + exp(-aLpha:*timeVari_t) - exp(-aLpha:*abs(timeVari_t-timeVari_s))) /// + 5*aLpha:*(exp(-aLpha:*timeVari_s):*timeVari_s + exp(-aLpha:*timeVari_t):*timeVari_t - exp(-aLpha:*abs(timeVari_t-timeVari_s)):*abs(timeVari_t-timeVari_s)) /// + aLpha^2:*(exp(-aLpha:*timeVari_s):*timeVari_s:*timeVari_s + exp(-aLpha:*timeVari_t):*timeVari_t:*timeVari_t - exp(-aLpha:*abs(timeVari_t-timeVari_s)):*abs(timeVari_t-timeVari_s):*abs(timeVari_t-timeVari_s)) /// + 8*aLpha:*minMatrix) )) // DIFFERENTIATE WITH RESPECT TO eta AND tau p_Rdotdot[2] = &(makesymmetric( /// (tau/aLpha^3):*(3:*(oneMatrix - exp(-aLpha:*timeVari_s) - exp(-aLpha:*timeVari_t) + exp(-aLpha:*abs(timeVari_t-timeVari_s))) /// + aLpha:*(-exp(-aLpha:*timeVari_s):*timeVari_s - exp(-aLpha:*timeVari_t):*timeVari_t + exp(-aLpha:*abs(timeVari_t-timeVari_s)):*abs(timeVari_t-timeVari_s) - 4:*minMatrix)) )) // DIFFERENTIATE WITH RESPECT TO tau AND tau p_Rdotdot[3] = &(makesymmetric( /// (1/aLpha^3):*(2*aLpha:*minMatrix + exp(-aLpha:*timeVari_s) + exp(-aLpha:*timeVari_t) - oneMatrix - exp(-aLpha:*abs(timeVari_t-timeVari_s))) )) } // END OF FUNCTION Rdotdot_et() void Rdot_eo(real scalar ni, real rowvector Rparameters, /// real matrix timeVari_s, real matrix timeVari_t, real matrix minMatrix, /// pointer(real matrix) colvector p_Rdot) { real scalar aLpha, omega real matrix oneMatrix aLpha = exp(Rparameters[1]) omega = Rparameters[2] oneMatrix = J(ni,ni,1) p_Rdot = J(2,1,NULL) // DIFFERENTIATE WITH RESPECT TO eta p_Rdot[1] = &(makesymmetric( /// (omega/(2*aLpha)):*(oneMatrix + exp(-aLpha:*abs(timeVari_t-timeVari_s)):*(oneMatrix + aLpha:*abs(timeVari_t-timeVari_s)) /// - exp(-aLpha:*timeVari_s):*(oneMatrix + aLpha:*timeVari_s) - exp(-aLpha:*timeVari_t):*(oneMatrix + aLpha:*timeVari_t)) )) // DIFFERENTIATE WITH RESPECT TO omega p_Rdot[2] = &(makesymmetric( /// (1/(2*aLpha)):*(2*aLpha:*minMatrix + exp(-aLpha:*timeVari_s) + exp(-aLpha:*timeVari_t) - oneMatrix - exp(-aLpha:*abs(timeVari_t-timeVari_s))) )) } // END OF FUNCTION Rdot_eo() void Rdotdot_eo(real scalar ni, real rowvector Rparameters, /// real matrix timeVari_s, real matrix timeVari_t, real matrix minMatrix, /// pointer(real matrix) colvector p_Rdotdot) { real scalar aLpha, omega real matrix oneMatrix aLpha = exp(Rparameters[1]) omega = Rparameters[2] oneMatrix = J(ni,ni,1) p_Rdotdot = J(3,1,NULL) // DIFFERENTIATE WITH RESPECT TO eta AND eta p_Rdotdot[1] = &(makesymmetric( /// (-omega/(2*aLpha)):*exp(-aLpha:*(timeVari_s+timeVari_t+abs(timeVari_t-timeVari_s))):*(exp(aLpha:*(timeVari_s+timeVari_t+abs(timeVari_t-timeVari_s))) /// + exp(aLpha:*(timeVari_s+timeVari_t)):*(oneMatrix + aLpha:*abs(timeVari_t-timeVari_s) + aLpha^2:*abs(timeVari_t-timeVari_s):*abs(timeVari_t-timeVari_s)) /// - exp(aLpha:*(timeVari_t+abs(timeVari_t-timeVari_s))):*(oneMatrix + aLpha:*timeVari_s + aLpha^2:*timeVari_s:*timeVari_s) /// - exp(aLpha:*(timeVari_s+abs(timeVari_t-timeVari_s))):*(oneMatrix + aLpha:*timeVari_t + aLpha^2:*timeVari_t:*timeVari_t)) )) // DIFFERENTIATE WITH RESPECT TO eta AND omega p_Rdotdot[2] = &(makesymmetric( /// (1/(2*aLpha)):*(oneMatrix + exp(-aLpha:*abs(timeVari_t-timeVari_s)):*(oneMatrix + aLpha:*abs(timeVari_t-timeVari_s)) /// - exp(-aLpha:*timeVari_s):*(oneMatrix + aLpha:*timeVari_s) - exp(-aLpha:*timeVari_t):*(oneMatrix+aLpha:*timeVari_t)) )) // DIFFERENTIATE WITH RESPECT TO omega AND omega p_Rdotdot[3] = &(J(ni,ni,0)) } // END OF FUNCTION Rdotdot_eo() void Rdot_it(real scalar ni, real rowvector Rparameters, /// real matrix timeVari_s, real matrix timeVari_t, real matrix minMatrix, /// pointer(real matrix) colvector p_Rdot) { real scalar aLpha, tau real matrix oneMatrix aLpha = sqrt(1/Rparameters[1]) tau = Rparameters[2] oneMatrix = J(ni,ni,1) p_Rdot = J(2,1,NULL) // DIFFERENTIATE WITH RESPECT TO iota p_Rdot[1] = &(makesymmetric( /// 0.25*(tau^2):*(exp(-aLpha:*timeVari_s):*timeVari_s + exp(-aLpha:*timeVari_t):*timeVari_t - (3/aLpha):*(oneMatrix - exp(-aLpha:*timeVari_s) - exp(-aLpha:*timeVari_t) + exp(-aLpha:*abs(timeVari_t-timeVari_s))) /// - exp(-aLpha:*abs(timeVari_t-timeVari_s)):*abs(timeVari_t-timeVari_s) + 4:*minMatrix) )) // DIFFERENTIATE WITH RESPECT TO tau p_Rdot[2] = &(makesymmetric( /// tau/(aLpha^3):*(2*aLpha:*minMatrix + exp(-aLpha:*timeVari_s) + exp(-aLpha:*timeVari_t) - oneMatrix - exp(-aLpha:*abs(timeVari_t-timeVari_s))) )) } // END OF FUNCTION Rdot_it() void Rdotdot_it(real scalar ni, real rowvector Rparameters, /// real matrix timeVari_s, real matrix timeVari_t, real matrix minMatrix, /// pointer(real matrix) colvector p_Rdotdot) { real scalar aLpha, tau real matrix oneMatrix aLpha = sqrt(1/Rparameters[1]) tau = Rparameters[2] oneMatrix = J(ni,ni,1) p_Rdotdot = J(3,1,NULL) // DIFFERENTIATE WITH RESPECT TO iota AND iota p_Rdotdot[1] = &(makesymmetric( /// -(((tau^2)*(aLpha^3))/8):*exp(-aLpha:*(timeVari_s+timeVari_t+abs(timeVari_t-timeVari_s))):* /// (- exp(aLpha:*(timeVari_t+abs(timeVari_t-timeVari_s))):*(timeVari_s:*timeVari_s+(3/aLpha^2):*oneMatrix+(3/aLpha):*timeVari_s) /// - exp(aLpha:*(timeVari_s+abs(timeVari_t-timeVari_s))):*(timeVari_t:*timeVari_t+(3/aLpha^2):*oneMatrix+(3/aLpha):*timeVari_t) /// + (3/aLpha^2):*exp(aLpha:*(timeVari_s+timeVari_t)) + (3/aLpha^2):*exp(aLpha:*(timeVari_s+timeVari_t+abs(timeVari_t-timeVari_s))) /// + (3/aLpha):*exp(aLpha:*(timeVari_s+timeVari_t)):*abs(timeVari_t-timeVari_s) + exp(aLpha:*(timeVari_s+timeVari_t)):*abs(timeVari_t-timeVari_s):*abs(timeVari_t-timeVari_s)) )) // DIFFERENTIATE WITH RESPECT TO iota AND tau p_Rdotdot[2] = &(makesymmetric( /// 0.5*tau:*(exp(-aLpha:*timeVari_s):*timeVari_s + exp(-aLpha:*timeVari_t):*timeVari_t - (3/aLpha):*(oneMatrix - exp(-aLpha:*timeVari_s) - exp(-aLpha:*timeVari_t) + exp(-aLpha:*abs(timeVari_t-timeVari_s))) /// - exp(-aLpha:*abs(timeVari_t-timeVari_s)):*abs(timeVari_t-timeVari_s) + 4:*minMatrix) )) // DIFFERENTIATE WITH RESPECT TO tau AND tau p_Rdotdot[3] = &(makesymmetric( /// 1/(aLpha^3):*(2*aLpha:*minMatrix + exp(-aLpha:*timeVari_s) + exp(-aLpha:*timeVari_t) - oneMatrix - exp(-aLpha:*abs(timeVari_t-timeVari_s))) )) } // END OF FUNCTION Rdotdot_it() void Rdot_io(real scalar ni, real rowvector Rparameters, /// real matrix timeVari_s, real matrix timeVari_t, real matrix minMatrix, /// pointer(real matrix) colvector p_Rdot) { real scalar aLpha, omega real matrix oneMatrix aLpha = sqrt(1/Rparameters[1]) omega = Rparameters[2] oneMatrix = J(ni,ni,1) p_Rdot = J(2,1,NULL) // DIFFERENTIATE WITH RESPECT TO iota p_Rdot[1] = &(makesymmetric( /// -0.25*aLpha*omega:*exp(-aLpha:*(timeVari_s+timeVari_t+abs(timeVari_t-timeVari_s))):*(exp(aLpha:*(timeVari_s+timeVari_t)):*(oneMatrix+aLpha:*abs(timeVari_t-timeVari_s)) /// + exp(aLpha:*(timeVari_s+timeVari_t+abs(timeVari_t-timeVari_s))) - exp(aLpha:*(timeVari_t+abs(timeVari_t-timeVari_s))):*(oneMatrix+aLpha:*timeVari_s) /// - exp(aLpha:*(timeVari_s+abs(timeVari_t-timeVari_s))):*(oneMatrix+aLpha:*timeVari_t)) )) // DIFFERENTIATE WITH RESPECT TO omega p_Rdot[2] = &(makesymmetric( /// (1/(2*aLpha)):*(2*aLpha:*minMatrix + exp(-aLpha:*timeVari_s) + exp(-aLpha:*timeVari_t) - oneMatrix - exp(-aLpha:*abs(timeVari_t-timeVari_s))) )) } // END OF FUNCTION Rdot_io() void Rdotdot_io(real scalar ni, real rowvector Rparameters, /// real matrix timeVari_s, real matrix timeVari_t, real matrix minMatrix, /// pointer(real matrix) colvector p_Rdotdot) { real scalar aLpha, omega real matrix oneMatrix aLpha = sqrt(1/Rparameters[1]) omega = Rparameters[2] oneMatrix = J(ni,ni,1) p_Rdotdot = J(3,1,NULL) // DIFFERENTIATE WITH RESPECT TO iota AND iota p_Rdotdot[1] = &(makesymmetric( /// -(omega/8)*(aLpha^5):*exp(-aLpha:*(timeVari_s+timeVari_t+abs(timeVari_t-timeVari_s))):*(exp(aLpha:*(timeVari_s+timeVari_t)):*abs(timeVari_t-timeVari_s):*abs(timeVari_t-timeVari_s) + /// exp(aLpha:*(timeVari_t+abs(timeVari_t-timeVari_s))):*(-timeVari_s:*timeVari_s+(1/aLpha^2):*oneMatrix+(1/aLpha):*timeVari_s) + /// exp(aLpha:*(timeVari_s+abs(timeVari_t-timeVari_s))):*(-timeVari_t:*timeVari_t+(1/aLpha^2):*oneMatrix+(1/aLpha):*timeVari_t) - /// (1/aLpha^2):*(exp(aLpha:*(timeVari_s+timeVari_t))+exp(aLpha:*(timeVari_s+timeVari_t+abs(timeVari_t-timeVari_s)))) - /// (1/aLpha):*exp(aLpha:*(timeVari_s+timeVari_t)):*abs(timeVari_t-timeVari_s)) )) // DIFFERENTIATE WITH RESPECT TO iota AND omega p_Rdotdot[2] = &(makesymmetric( /// -0.25*aLpha:*exp(-aLpha:*(timeVari_s+timeVari_t+abs(timeVari_t-timeVari_s))):*(exp(aLpha:*(timeVari_s+timeVari_t)):*(oneMatrix+aLpha:*abs(timeVari_t-timeVari_s)) /// + exp(aLpha:*(timeVari_s+timeVari_t+abs(timeVari_t-timeVari_s))) - exp(aLpha:*(timeVari_t+abs(timeVari_t-timeVari_s))):*(oneMatrix+aLpha:*timeVari_s) /// - exp(aLpha:*(timeVari_s+abs(timeVari_t-timeVari_s))):*(oneMatrix+aLpha:*timeVari_t)) )) // DIFFERENTIATE WITH RESPECT TO omega AND omega p_Rdotdot[3] = &(J(ni,ni,0)) } // END OF FUNCTION Rdotdot_io() void Rdot_ek(real scalar ni, real rowvector Rparameters, /// real matrix timeVari_s, real matrix timeVari_t, real matrix minMatrix, /// pointer(real matrix) colvector p_Rdot) { real scalar aLpha, tau, sigma real matrix oneMatrix aLpha = exp(Rparameters[1]) tau = exp(Rparameters[2]) sigma = exp(Rparameters[3]) oneMatrix = J(ni,ni,1) p_Rdot = J(3,1,NULL) // DIFFERENTIATE WITH RESPECT TO eta p_Rdot[1] = &(makesymmetric( /// (tau^2/(2*(aLpha^3))):*(3:*(oneMatrix - exp(-aLpha:*timeVari_s) - exp(-aLpha:*timeVari_t) + exp(-aLpha:*abs(timeVari_t-timeVari_s))) /// + aLpha:*(-exp(-aLpha:*timeVari_s):*timeVari_s - exp(-aLpha:*timeVari_t):*timeVari_t + exp(-aLpha:*abs(timeVari_t-timeVari_s)):*abs(timeVari_t-timeVari_s) - 4:*minMatrix)) )) // DIFFERENTIATE WITH RESPECT TO kappa = lntau p_Rdot[2] = &(makesymmetric( /// (tau^2/aLpha^3):*(2*aLpha:*minMatrix + exp(-aLpha:*timeVari_s) + exp(-aLpha:*timeVari_t) - oneMatrix - exp(-aLpha:*abs(timeVari_t-timeVari_s))) )) // DIFFERENTIATE WITH RESPECT TO lnsigma p_Rdot[3] = &(makesymmetric( /// 2*(sigma^2):*I(ni) )) } // END OF FUNCTION Rdot_ek() void Rdotdot_ek(real scalar ni, real rowvector Rparameters, /// real matrix timeVari_s, real matrix timeVari_t, real matrix minMatrix, /// pointer(real matrix) colvector p_Rdotdot) { real scalar aLpha, tau, sigma real matrix oneMatrix aLpha = exp(Rparameters[1]) tau = exp(Rparameters[2]) sigma = exp(Rparameters[3]) oneMatrix = J(ni,ni,1) p_Rdotdot = J(6,1,NULL) // DIFFERENTIATE WITH RESPECT TO eta AND eta p_Rdotdot[1] = &(makesymmetric( /// (tau^2/(2*(aLpha)^3)):*(9:*(-oneMatrix + exp(-aLpha:*timeVari_s) + exp(-aLpha:*timeVari_t) - exp(-aLpha:*abs(timeVari_t-timeVari_s))) /// + 5*aLpha:*(exp(-aLpha:*timeVari_s):*timeVari_s + exp(-aLpha:*timeVari_t):*timeVari_t - exp(-aLpha:*abs(timeVari_t-timeVari_s)):*abs(timeVari_t-timeVari_s)) /// + aLpha^2:*(exp(-aLpha:*timeVari_s):*timeVari_s:*timeVari_s + exp(-aLpha:*timeVari_t):*timeVari_t:*timeVari_t - exp(-aLpha:*abs(timeVari_t-timeVari_s)):*abs(timeVari_t-timeVari_s):*abs(timeVari_t-timeVari_s)) /// + 8*aLpha:*minMatrix) )) // DIFFERENTIATE WITH RESPECT TO eta AND kappa p_Rdotdot[2] = &(makesymmetric( /// (tau^2/aLpha^3):*(3:*(oneMatrix - exp(-aLpha:*timeVari_s) - exp(-aLpha:*timeVari_t) + exp(-aLpha:*abs(timeVari_t-timeVari_s))) /// + aLpha:*(-exp(-aLpha:*timeVari_s):*timeVari_s - exp(-aLpha:*timeVari_t):*timeVari_t + exp(-aLpha:*abs(timeVari_t-timeVari_s)):*abs(timeVari_t-timeVari_s) - 4:*minMatrix)) )) // DIFFERENTIATE WITH RESPECT TO eta AND lnsigma p_Rdotdot[3] = &(J(ni,ni,0)) // DIFFERENTIATE WITH RESPECT TO kappa AND kappa p_Rdotdot[4] = &(makesymmetric( /// (2*tau^2/aLpha^3):*(2*aLpha:*minMatrix + exp(-aLpha:*timeVari_s) + exp(-aLpha:*timeVari_t) - oneMatrix - exp(-aLpha:*abs(timeVari_t-timeVari_s))) )) // DIFFERENTIATE WITH RESPECT TO kappa AND lnsigma p_Rdotdot[5] = &(J(ni,ni,0)) // DIFFERENTIATE WITH RESPECT TO lnsigma AND lnsigma p_Rdotdot[6] = &(makesymmetric( /// 4*(sigma^2):*I(ni) )) } // END OF FUNCTION Rdotdot_ek() void Rdot_phi(real scalar ni, real rowvector Rparameters, /// real matrix timeVari_s, real matrix timeVari_t, real matrix minMatrix, /// pointer(real matrix) colvector p_Rdot) { p_Rdot = J(1,1,NULL) // DIFFERENTIATE WITH RESPECT TO phi p_Rdot[1] = &(minMatrix) } // END OF FUNCTION Rdot_phi() void Rdotdot_phi(real scalar ni, real rowvector Rparameters, /// real matrix timeVari_s, real matrix timeVari_t, real matrix minMatrix, /// pointer(real matrix) colvector p_Rdotdot) { p_Rdotdot = J(1,1,NULL) // DIFFERENTIATE WITH RESPECT TO phi AND phi p_Rdotdot[1] = &(J(ni,ni,0)) } // END OF FUNCTION Rdotdot_phi() void Rdot_lnphi(real scalar ni, real rowvector Rparameters, /// real matrix timeVari_s, real matrix timeVari_t, real matrix minMatrix, /// pointer(real matrix) colvector p_Rdot) { real scalar phi, sigma real matrix oneMatrix phi = exp(Rparameters[1]) sigma = exp(Rparameters[2]) oneMatrix = J(ni,ni,1) p_Rdot = J(2,1,NULL) // DIFFERENTIATE WITH RESPECT TO ln(phi) p_Rdot[1] = &(makesymmetric(phi:*minMatrix)) // DIFFERENTIATE WITH RESPECT TO lnsigma p_Rdot[2] = &(makesymmetric( /// 2*(sigma^2):*I(ni) )) } // END OF FUNCTION Rdot_lnphi() void Rdotdot_lnphi(real scalar ni, real rowvector Rparameters, /// real matrix timeVari_s, real matrix timeVari_t, real matrix minMatrix, /// pointer(real matrix) colvector p_Rdotdot) { real scalar phi, tau, sigma real matrix oneMatrix phi = exp(Rparameters[1]) sigma = exp(Rparameters[2]) oneMatrix = J(ni,ni,1) p_Rdotdot = J(3,1,NULL) // DIFFERENTIATE WITH RESPECT TO lnphi AND lnphi p_Rdotdot[1] = &(makesymmetric(phi:*minMatrix)) // DIFFERENTIATE WITH RESPECT TO lnphi AND lnsigma p_Rdotdot[2] = &(J(ni,ni,0)) // DIFFERENTIATE WITH RESPECT TO lnsigma AND lnsigma p_Rdotdot[3] = &(makesymmetric( /// 4*(sigma^2):*I(ni) )) } // END OF FUNCTION Rdotdot_lnphi() void Rconversion(real scalar inputRstyle, real scalar inputIOU1, real scalar inputIOU2, /// real scalar outputRstyle, real scalar outputIOU1, real scalar outputIOU2) { real scalar aLpha, tau // CONVERT inputIOU1 -> aLpha if(inputRstyle==3 | inputRstyle==4 | inputRstyle==7) { // eta -> alpha aLpha = exp(inputIOU1) } else if (inputRstyle==5 | inputRstyle==6) { // iota -> alpha aLpha = sqrt(1/inputIOU1) } else { aLpha = inputIOU1 } // CONVERT inputIOU2 -> tau if (inputRstyle==2 | inputRstyle==4 | inputRstyle==6) { // omega -> tau tau = sqrt((aLpha^2)*inputIOU2) } else if(inputRstyle==7) { // kappa -> tau tau = exp(inputIOU2) } else { tau = inputIOU2 } // CONVERT aLpha -> outputIOU1 if (outputRstyle==3 | outputRstyle==4 | outputRstyle==7) { // alpha -> eta outputIOU1 = ln(aLpha) } else if (outputRstyle==5 | outputRstyle==6) { // aLpha -> iota outputIOU1 = 1/(aLpha^2) } else { outputIOU1 = aLpha } // CONVERT tau -> outputIOU2 if (outputRstyle==2 | outputRstyle==4 | outputRstyle==6) { // tau -> omega outputIOU2 = (tau/aLpha)^2 } else if (outputRstyle==7) { // tau -> kappa outputIOU2 = ln(tau) } else { outputIOU2 = tau } } // END OF FUNCTION Rconversion() // SEE "xtIOU/do files/tests delta method for sigma and sigmaSquared.do" FOR CHECKING OF MY IMPLEMENTATION OF THE DELTA METHOD void tabTheta_sigmaSq(real rowvector transParameters, real matrix transVariance, /// real scalar colsZ, real scalar numberOfGparas, /// real matrix table_theta) { real scalar numberOfParameters, fnpara, fnrow, fncol, dpara, drow, dcol real matrix variance, fnderiv, positions numberOfParameters = cols(transParameters) positions = invvech(range(1, numberOfGparas, 1)) table_theta = J(numberOfParameters,4,.) table_theta[|1,2 \ numberOfParameters,2|] = sqrt(diagonal(transVariance)) // Std Err FOR TRANSFORMED PARAMETERS fnderiv = J(numberOfParameters,numberOfParameters,0) fncol = 1 fnrow = 0 for(fnpara=1; fnpara<=numberOfGparas; fnpara++) { fnrow = fnrow + 1 if(fnrow > colsZ) { fncol = fncol + 1 fnrow = fncol } // CALCULATES CONFIDENCE INTERVALS FOR VARIANCE PARAMETERS (i.e. DIAGONAL ENTRIES) OF MATRIX G if(fnrow==fncol) { table_theta[fnpara,1] = exp(2*transParameters[fnpara]) table_theta[fnpara,3] = exp(2*(transParameters[fnpara] - 1.96*table_theta[fnpara,2])) table_theta[fnpara,4] = exp(2*(transParameters[fnpara] + 1.96*table_theta[fnpara,2])) } else { table_theta[fnpara,1] = tanh(transParameters[fnpara])*exp(transParameters[positions[fnrow,fnrow]])*exp(transParameters[positions[fncol,fncol]]) } // CALCULATES DERIVATIVE OF TRANSFORMATION FUNCTION (transParameters -> theta) FOR THE PARAMETERS OF MATRIX G dcol = 1 drow = 0 for(dpara=1; dpara<=numberOfGparas; dpara++) { drow = drow + 1 if(drow > colsZ) { dcol = dcol + 1 drow = dcol } if(fnpara==dpara) { if(fnrow==fncol) { fnderiv[fnpara,fnpara] = 2*exp(2*transParameters[fnpara]) } else { fnderiv[fnpara,fnpara] = (1/cosh(transParameters[fnpara]))^2*exp(transParameters[positions[fnrow,fnrow]])*exp(transParameters[positions[fncol,fncol]]) } } else { if((fncol==drow & drow==dcol) | (fnrow==drow & drow==dcol)) { fnderiv[fnpara,dpara] = tanh(transParameters[fnpara])*exp(transParameters[positions[fnrow,fnrow]])*exp(transParameters[positions[fncol,fncol]]) } // ELSE A ZERO VALUE } } // END OF dpara FOR-LOOP } // END OF fnpara FOR-LOOP // CALCULATES DERIVATIVE OF TRANSFORMATION FUNCTION FOR THE PARAMETERS OF MATRIX W // CALCULATES STANDARD ERRORS AND CONFIDENCE INTERVALS FOR THE PARAMETERS OF MATRIX W for(fnpara=numberOfGparas+1; fnpara<=numberOfParameters-1; fnpara++) { fnderiv[fnpara,fnpara] = exp(transParameters[fnpara]) // e.g. DERIVATIVE WHEN CONVERTING FROM lnalpha TO aLpha table_theta[fnpara,1] = exp(transParameters[fnpara]) table_theta[fnpara,3] = exp(transParameters[fnpara] - 1.96*table_theta[fnpara,2]) table_theta[fnpara,4] = exp(transParameters[fnpara] + 1.96*table_theta[fnpara,2]) } // CALCULATE DERIVATIVE OF TRANSFORMATION FUNCTION FOR lnsigma -> sigmaSquared // CALCULATES STANDARD ERRORS AND CONFIDENCE INTERVALS FOR sigmaSquared fnderiv[numberOfParameters,numberOfParameters] = 2*exp(2*transParameters[numberOfParameters]) table_theta[numberOfParameters,1] = exp(2*transParameters[numberOfParameters]) table_theta[numberOfParameters,3] = exp(2*(transParameters[numberOfParameters] - 1.96*table_theta[numberOfParameters,2])) table_theta[numberOfParameters,4] = exp(2*(transParameters[numberOfParameters] + 1.96*table_theta[numberOfParameters,2])) // COVARIANCE MATRIX OF UNTRANSFORMED PARAMETERS variance = fnderiv*transVariance*fnderiv' table_theta[|1,2 \ numberOfParameters,2|] = sqrt(diagonal(variance)) // Std Err FOR UNTRANSFORMED PARAMETERS // CALCULATES CONFIDENCE INTERVALS FOR THE COVARIANCES (i.e. OFF-DIAGONAL ELEMENTS) OF MATRIX G for(fnpara=1; fnpara<=numberOfGparas; fnpara++) { if(table_theta[fnpara,3]==.) { table_theta[fnpara,3] = table_theta[fnpara,1] - 1.96*table_theta[fnpara,2] table_theta[fnpara,4] = table_theta[fnpara,1] + 1.96*table_theta[fnpara,2] } } } // END OF tabTheta_sigmaSq // INPUT: ESTIMATES WITH CORRESPONDING COVARIANCE MATRIX // IF ZSCORES AND PVALUES ARE REQUIRED THEN SET style TO 1 // OUTPUT: MATRIX CONTAINING ESTIMATES, SEs, [ZSCORES, PVALUES], LCIs, UCIs void tabulatingBeta(real colvector vec_estimates, real matrix mat_covariance, real matrix table_beta) { real colvector blank blank = J(rows(vec_estimates),1,0) // COEFFICIENTS AND STANDARD ERRORS table_beta = vec_estimates table_beta = table_beta, sqrt(diagonal(mat_covariance)) // Z-SCORES AND P-VALUES table_beta = table_beta, blank table_beta[,3] = table_beta[,1] :/ table_beta[,2] table_beta = table_beta, blank table_beta[,4] = 2*(1 :- normal(abs(table_beta[,3]))) // ADD lci table_beta = table_beta, blank table_beta[,5] = table_beta[,1] :- invnormal(0.975):*table_beta[,2] // ADD uci table_beta = table_beta, blank table_beta[,6] = table_beta[,1] :+ invnormal(0.975):*table_beta[,2] } // END OF FUNCTION tabulating() end