R/gradient_stochastic_adjustment.R

In diseq: Estimation Methods for Markets in Equilibrium and Disequilibrium

#' @include system_stochastic_adjustment.R

setMethod("calculate_system_scores",
          signature(object = "system_stochastic_adjustment"),
          function(object) {
  # nolint start
  sd <- object@demand@sigma
  ss <- object@supply@sigma
  vd <- object@demand@var
  vs <- object@supply@var
  ad <- object@demand@alpha
  as <- object@supply@alpha
  bd <- object@demand@beta
  bs <- object@supply@beta
  xd <- object@demand@control_matrix
  xs <- object@supply@control_matrix
  r <- object@rho
  rds <- object@rho_ds
  rdp <- object@rho_dp
  rsp <- object@rho_sp
  gm <- object@gamma
  dl <- object@delta
  LP <- object@lagged_price_vector
  xp <- object@price_equation@control_matrix
  bp <- object@price_equation@beta
  sp <- object@price_equation@sigma
  vp <- object@price_equation@var
  psiD <- object@psi_D
  psiS <- object@psi_S
  PsiD <- object@Psi_D
  PsiS <- object@Psi_S
  rDP <- object@rho_DP
  rDS <- object@rho_DS
  rSP <- object@rho_SP
  zeta <- object@zeta
  zetaDS <- object@zeta_DS
  zetaDP <- object@zeta_DP
  zetaSP <- object@zeta_SP
  zetaDD <- object@zeta_DD
  zetaSS <- object@zeta_SS
  sD <- object@sigma_D
  sS <- object@sigma_S
  sP <- object@sigma_P
  vD <- object@var_D
  vS <- object@var_S
  vP <- object@var_P
  hD <- object@h_D
  hS <- object@h_S
  hP <- object@h_P
  zDP <- object@z_DP
  zPD <- object@z_PD
  zSP <- object@z_SP
  zPS <- object@z_PS
  omegaD <- object@omega_D
  omegaS <- object@omega_S
  wD <- object@w_D
  wS <- object@w_S
  gD <- object@g_D
  gS <- object@g_S
  LD <- object@L_D
  LS <- object@L_S
  xbd <- xd %*% bd
  xbs <- xs %*% bs
  xbp <- xp %*% bp

  palpha_d <- -(LD*zetaDD**(15/2)*(dl*sD**2*sP*sS**2*vD*vP*vS*zeta**2*zetaSS**7*(gS*zetaSP - zPD*zeta)*(LP*gm + xbd + xbp*gm - xbs) + dl*sD*vP*sS**2*vD*vP*vS*zeta**2*zetaSS**7*(ad + dl)*(gS*zetaDS - zDP*zeta)*(LP*gm + xbd + xbp*gm - xbs) - gS*sD*sP*zetaSS**(15/2)*(as*dl*sD*sP*sS*vD*vP*vS*zeta**2*sqrt(zetaSS)*(LP*gm + xbd + xbp*gm - xbs) + as*hS*sD*sP*vD*vP*vS*zeta**2*sqrt(zetaSS)*(as*(gm**2*vp + 2*gm*sp*(rdp*sd - rsp*ss) - 2*rds*sd*ss + vd + vs) + dl*ss*(gm*rsp*sp + rds*sd - ss)) + sD*sS*vD*(omegaS*zetaSP - zPD*zeta**2)*(rSP*sP*sS*(as*vP*(as*(gm**2*vp + 2*gm*sp*(rdp*sd - rsp*ss) - 2*rds*sd*ss + vd + vs) + dl*ss*(gm*rsp*sp + rds*sd - ss)) + vS*(gm**2*vp + 2*gm*sp*(rdp*sd - rsp*ss) - 2*rds*sd*ss + vd + vs)) - vP*vS*(2*as*(gm**2*vp + 2*gm*sp*(rdp*sd - rsp*ss) - 2*rds*sd*ss + vd + vs) + dl*ss*(gm*rsp*sp + rds*sd - ss))) + sP*sS*vP*(omegaS*zetaDS - zDP*zeta**2)*(rDS*sD*sS*(as*vD*(as*(gm**2*vp + 2*gm*sp*(rdp*sd - rsp*ss) - 2*rds*sd*ss + vd + vs) + dl*ss*(gm*rsp*sp + rds*sd - ss)) + vS*(ad**2*(gm**2*vp + 2*gm*sp*(rdp*sd - rsp*ss) - 2*rds*sd*ss + vd + vs) + ad*dl*(gm**2*vp + 2*gm*sp*(rdp*sd - rsp*ss) - 2*rds*sd*ss + vd + sd*(gm*rdp*sp - rds*ss + sd) + vs) + dl**2*sd*(gm*rdp*sp - rds*ss + sd))) - vD*vS*(2*ad*as*(gm**2*vp + 2*gm*sp*(rdp*sd - rsp*ss) - 2*rds*sd*ss + vd + vs) + dl**2*ss*(gm*rsp*sp + rds*sd - ss) + dl*(ad*ss*(gm*rsp*sp + rds*sd - ss) + as*sd*(gm*rdp*sp - rds*ss + sd) + as*(gm**2*vp + 2*gm*sp*(rdp*sd - rsp*ss) - 2*rds*sd*ss + vd + vs))))) + sD**2*sS**2*vD*vP*vS*zeta**2*zetaSS**(11/2)*(gS*hP*zetaSP*zetaSS**(3/2) - hD*zDP*zeta*zetaSS**(3/2) + zeta*zetaSS**2 + zeta*zetaSS*(hD**2 - hP**2))*(gm**2*vp + 2*gm*sp*(rdp*sd - rsp*ss) - 2*rds*sd*ss + vd + vs) - sD*sP*sS**2*vS*zetaSS**5*(gS*zetaSS*(omegaS*rDP*zeta**2*sqrt(zetaSS) + omegaS*zetaDP*zetaSS**(3/2) + zeta**2*zetaSS*(hD*rSP + hP*rDS - 2*hS*rDP)) - zeta**3*zetaSS**(3/2)*(hD*hP + 2*rDP*wD + rDP))*(rDP*sD*sP*(vD*(gm**2*vp + 2*gm*sp*(rdp*sd - rsp*ss) - 2*rds*sd*ss + vd + vs) + vP*(ad**2*(gm**2*vp + 2*gm*sp*(rdp*sd - rsp*ss) - 2*rds*sd*ss + vd + vs) + ad*dl*(gm**2*vp + 2*gm*sp*(rdp*sd - rsp*ss) - 2*rds*sd*ss + vd + sd*(gm*rdp*sp - rds*ss + sd) + vs) + dl**2*sd*(gm*rdp*sp - rds*ss + sd))) - vD*vP*(2*ad*(gm**2*vp + 2*gm*sp*(rdp*sd - rsp*ss) - 2*rds*sd*ss + vd + vs) + dl*(gm**2*vp + gm*rdp*sd*sp + 2*gm*sp*(rdp*sd - rsp*ss) - 3*rds*sd*ss + 2*vd + vs))) + vP*sS**2*vD*vP*vS*zeta**2*zetaSS**(13/2)*(ad**2*(gm**2*vp + 2*gm*sp*(rdp*sd - rsp*ss) - 2*rds*sd*ss + vd + vs) + ad*dl*(gm**2*vp + 2*gm*sp*(rdp*sd - rsp*ss) - 2*rds*sd*ss + vd + sd*(gm*rdp*sp - rds*ss + sd) + vs) + dl**2*sd*(gm*rdp*sp - rds*ss + sd))*(gS*hD*zetaDS*sqrt(zetaSS) - hD*zDP*zeta*sqrt(zetaSS) + zeta*zetaSS)) + LS*zetaSS**(15/2)*(as*dl*sD**2*vP*sS*vD*vP*vS*zeta**2*zetaDD**7*(gD*zetaDS - zSP*zeta)*(LP*gm + xbd + xbp*gm - xbs) + as*sD**2*vP*vD*vP*vS*zeta**2*zetaDD**(13/2)*(as*(gm**2*vp + 2*gm*sp*(rdp*sd - rsp*ss) - 2*rds*sd*ss + vd + vs) + dl*ss*(gm*rsp*sp + rds*sd - ss))*(gD*hS*sqrt(zetaDD)*zetaDS - hS*zSP*zeta*sqrt(zetaDD) + zeta*zetaDD) + dl*sD**2*sP*sS**2*vD*vP*vS*zeta**2*zetaDD**7*(gD*zetaDP - zPS*zeta)*(LP*gm + xbd + xbp*gm - xbs) - gD*sP*sS*zetaDD**(15/2)*(dl*sD*sP*sS*vD*vP*vS*zeta**2*sqrt(zetaDD)*(ad + dl)*(LP*gm + xbd + xbp*gm - xbs) + hD*sP*sS*vD*vP*vS*zeta**2*sqrt(zetaDD)*(ad**2*(gm**2*vp + 2*gm*sp*(rdp*sd - rsp*ss) - 2*rds*sd*ss + vd + vs) + ad*dl*(gm**2*vp + 2*gm*sp*(rdp*sd - rsp*ss) - 2*rds*sd*ss + vd + sd*(gm*rdp*sp - rds*ss + sd) + vs) + dl**2*sd*(gm*rdp*sp - rds*ss + sd)) + sD*sP*vP*(omegaD*zetaDS - zSP*zeta**2)*(rDS*sD*sS*(as*vD*(as*(gm**2*vp + 2*gm*sp*(rdp*sd - rsp*ss) - 2*rds*sd*ss + vd + vs) + dl*ss*(gm*rsp*sp + rds*sd - ss)) + vS*(ad**2*(gm**2*vp + 2*gm*sp*(rdp*sd - rsp*ss) - 2*rds*sd*ss + vd + vs) + ad*dl*(gm**2*vp + 2*gm*sp*(rdp*sd - rsp*ss) - 2*rds*sd*ss + vd + sd*(gm*rdp*sp - rds*ss + sd) + vs) + dl**2*sd*(gm*rdp*sp - rds*ss + sd))) - vD*vS*(2*ad*as*(gm**2*vp + 2*gm*sp*(rdp*sd - rsp*ss) - 2*rds*sd*ss + vd + vs) + dl**2*ss*(gm*rsp*sp + rds*sd - ss) + dl*(ad*ss*(gm*rsp*sp + rds*sd - ss) + as*sd*(gm*rdp*sp - rds*ss + sd) + as*(gm**2*vp + 2*gm*sp*(rdp*sd - rsp*ss) - 2*rds*sd*ss + vd + vs)))) + sD*sS*vS*(omegaD*zetaDP - zPS*zeta**2)*(rDP*sD*sP*(vD*(gm**2*vp + 2*gm*sp*(rdp*sd - rsp*ss) - 2*rds*sd*ss + vd + vs) + vP*(ad**2*(gm**2*vp + 2*gm*sp*(rdp*sd - rsp*ss) - 2*rds*sd*ss + vd + vs) + ad*dl*(gm**2*vp + 2*gm*sp*(rdp*sd - rsp*ss) - 2*rds*sd*ss + vd + sd*(gm*rdp*sp - rds*ss + sd) + vs) + dl**2*sd*(gm*rdp*sp - rds*ss + sd))) - vD*vP*(2*ad*(gm**2*vp + 2*gm*sp*(rdp*sd - rsp*ss) - 2*rds*sd*ss + vd + vs) + dl*(gm**2*vp + gm*rdp*sd*sp + 2*gm*sp*(rdp*sd - rsp*ss) - 3*rds*sd*ss + 2*vd + vs)))) - sD**2*sP*sS*vD*zetaDD**5*(gD*zetaDD*(omegaD*rSP*zeta**2*sqrt(zetaDD) + omegaD*zetaDD**(3/2)*zetaSP + zeta**2*zetaDD*(-2*hD*rSP + hP*rDS + hS*rDP)) - zeta**3*zetaDD**(3/2)*(hP*hS + 2*rSP*wS + rSP))*(rSP*sP*sS*(as*vP*(as*(gm**2*vp + 2*gm*sp*(rdp*sd - rsp*ss) - 2*rds*sd*ss + vd + vs) + dl*ss*(gm*rsp*sp + rds*sd - ss)) + vS*(gm**2*vp + 2*gm*sp*(rdp*sd - rsp*ss) - 2*rds*sd*ss + vd + vs)) - vP*vS*(2*as*(gm**2*vp + 2*gm*sp*(rdp*sd - rsp*ss) - 2*rds*sd*ss + vd + vs) + dl*ss*(gm*rsp*sp + rds*sd - ss))) + sD**2*sS**2*vD*vP*vS*zeta**2*zetaDD**(11/2)*(gD*hP*zetaDD**(3/2)*zetaDP - hS*zSP*zeta*zetaDD**(3/2) + zeta*zetaDD**2 - zeta*zetaDD*(hP**2 - hS**2))*(gm**2*vp + 2*gm*sp*(rdp*sd - rsp*ss) - 2*rds*sd*ss + vd + vs)))/(dl**3*sD**2*vP*sS**2*vD*vP*vS*zeta**3*zetaDD**(15/2)*zetaSS**(15/2))
  pbeta_d1 <- sweep(-xd, MARGIN = 1, (LD*zetaDD**2*zetaSS**(3/2)*(-as*gS*sD*sP*zetaSS + sD*sS*(gS*zetaSP - zPD*zeta) + sP*sS*(ad + dl)*(gS*zetaDS - zDP*zeta)) + LS*zetaDD**(3/2)*zetaSS**2*(as*sD*sP*(gD*zetaDS - zSP*zeta) - gD*sP*sS*zetaDD*(ad + dl) + sD*sS*(gD*zetaDP - zPS*zeta)))/(dl*sD*sP*sS*zeta*zetaDD**2*zetaSS**2), `*`)
  palpha_s <- (LD*zetaDD**(15/2)*(ad*dl*sD*vP*sS**2*vD*vP*vS*zeta**2*zetaSS**7*(gS*zetaDS - zDP*zeta)*(LP*gm + xbd + xbp*gm - xbs) + ad*vP*sS**2*vD*vP*vS*zeta**2*zetaSS**(13/2)*(ad*(gm**2*vp + 2*gm*sp*(rdp*sd - rsp*ss) - 2*rds*sd*ss + vd + vs) + dl*sd*(gm*rdp*sp - rds*ss + sd))*(gS*hD*zetaDS*sqrt(zetaSS) - hD*zDP*zeta*sqrt(zetaSS) + zeta*zetaSS) + dl*sD**2*sP*sS**2*vD*vP*vS*zeta**2*zetaSS**7*(gS*zetaSP - zPD*zeta)*(LP*gm + xbd + xbp*gm - xbs) - gS*sD*sP*zetaSS**(15/2)*(dl*sD*sP*sS*vD*vP*vS*zeta**2*sqrt(zetaSS)*(as - dl)*(LP*gm + xbd + xbp*gm - xbs) - hS*sD*sP*vD*vP*vS*zeta**2*sqrt(zetaSS)*(as**2*(-gm**2*vp - 2*gm*sp*(rdp*sd - rsp*ss) + 2*rds*sd*ss - vd - vs) + as*dl*(gm**2*vp + 2*gm*sp*(rdp*sd - rsp*ss) - 2*rds*sd*ss + vd + vs - ss*(gm*rsp*sp + rds*sd - ss)) + dl**2*ss*(gm*rsp*sp + rds*sd - ss)) - sD*sS*vD*(omegaS*zetaSP - zPD*zeta**2)*(rSP*sP*sS*(vP*(as**2*(-gm**2*vp - 2*gm*sp*(rdp*sd - rsp*ss) + 2*rds*sd*ss - vd - vs) + as*dl*(gm**2*vp + 2*gm*sp*(rdp*sd - rsp*ss) - 2*rds*sd*ss + vd + vs - ss*(gm*rsp*sp + rds*sd - ss)) + dl**2*ss*(gm*rsp*sp + rds*sd - ss)) - vS*(gm**2*vp + 2*gm*sp*(rdp*sd - rsp*ss) - 2*rds*sd*ss + vd + vs)) + vP*vS*(2*as*(gm**2*vp + 2*gm*sp*(rdp*sd - rsp*ss) - 2*rds*sd*ss + vd + vs) - dl*(gm**2*vp - gm*rsp*sp*ss + 2*gm*sp*(rdp*sd - rsp*ss) - 3*rds*sd*ss + vd + 2*vs))) + sP*sS*vP*(omegaS*zetaDS - zDP*zeta**2)*(rDS*sD*sS*(ad*vS*(ad*(gm**2*vp + 2*gm*sp*(rdp*sd - rsp*ss) - 2*rds*sd*ss + vd + vs) + dl*sd*(gm*rdp*sp - rds*ss + sd)) - vD*(as**2*(-gm**2*vp - 2*gm*sp*(rdp*sd - rsp*ss) + 2*rds*sd*ss - vd - vs) + as*dl*(gm**2*vp + 2*gm*sp*(rdp*sd - rsp*ss) - 2*rds*sd*ss + vd + vs - ss*(gm*rsp*sp + rds*sd - ss)) + dl**2*ss*(gm*rsp*sp + rds*sd - ss))) - vD*vS*(2*ad*as*(gm**2*vp + 2*gm*sp*(rdp*sd - rsp*ss) - 2*rds*sd*ss + vd + vs) - dl**2*sd*(gm*rdp*sp - rds*ss + sd) + dl*(ad*ss*(gm*rsp*sp + rds*sd - ss) - ad*(gm**2*vp + 2*gm*sp*(rdp*sd - rsp*ss) - 2*rds*sd*ss + vd + vs) + as*sd*(gm*rdp*sp - rds*ss + sd))))) + sD**2*sS**2*vD*vP*vS*zeta**2*zetaSS**(11/2)*(gS*hP*zetaSP*zetaSS**(3/2) - hD*zDP*zeta*zetaSS**(3/2) + zeta*zetaSS**2 + zeta*zetaSS*(hD**2 - hP**2))*(gm**2*vp + 2*gm*sp*(rdp*sd - rsp*ss) - 2*rds*sd*ss + vd + vs) + sD*sP*sS**2*vS*zetaSS**5*(-gS*zetaSS*(omegaS*rDP*zeta**2*sqrt(zetaSS) + omegaS*zetaDP*zetaSS**(3/2) + zeta**2*zetaSS*(hD*rSP + hP*rDS - 2*hS*rDP)) + zeta**3*zetaSS**(3/2)*(hD*hP + 2*rDP*wD + rDP))*(rDP*sD*sP*(ad*vP*(ad*(gm**2*vp + 2*gm*sp*(rdp*sd - rsp*ss) - 2*rds*sd*ss + vd + vs) + dl*sd*(gm*rdp*sp - rds*ss + sd)) + vD*(gm**2*vp + 2*gm*sp*(rdp*sd - rsp*ss) - 2*rds*sd*ss + vd + vs)) - vD*vP*(2*ad*(gm**2*vp + 2*gm*sp*(rdp*sd - rsp*ss) - 2*rds*sd*ss + vd + vs) + dl*sd*(gm*rdp*sp - rds*ss + sd)))) - LS*zetaSS**(15/2)*(-dl*sD**2*vP*sS*vD*vP*vS*zeta**2*zetaDD**7*(as - dl)*(gD*zetaDS - zSP*zeta)*(LP*gm + xbd + xbp*gm - xbs) - dl*sD**2*sP*sS**2*vD*vP*vS*zeta**2*zetaDD**7*(gD*zetaDP - zPS*zeta)*(LP*gm + xbd + xbp*gm - xbs) + gD*sP*sS*zetaDD**(15/2)*(ad*dl*sD*sP*sS*vD*vP*vS*zeta**2*sqrt(zetaDD)*(LP*gm + xbd + xbp*gm - xbs) + ad*hD*sP*sS*vD*vP*vS*zeta**2*sqrt(zetaDD)*(ad*(gm**2*vp + 2*gm*sp*(rdp*sd - rsp*ss) - 2*rds*sd*ss + vd + vs) + dl*sd*(gm*rdp*sp - rds*ss + sd)) + sD*sP*vP*(omegaD*zetaDS - zSP*zeta**2)*(rDS*sD*sS*(ad*vS*(ad*(gm**2*vp + 2*gm*sp*(rdp*sd - rsp*ss) - 2*rds*sd*ss + vd + vs) + dl*sd*(gm*rdp*sp - rds*ss + sd)) - vD*(as**2*(-gm**2*vp - 2*gm*sp*(rdp*sd - rsp*ss) + 2*rds*sd*ss - vd - vs) + as*dl*(gm**2*vp + 2*gm*sp*(rdp*sd - rsp*ss) - 2*rds*sd*ss + vd + vs - ss*(gm*rsp*sp + rds*sd - ss)) + dl**2*ss*(gm*rsp*sp + rds*sd - ss))) - vD*vS*(2*ad*as*(gm**2*vp + 2*gm*sp*(rdp*sd - rsp*ss) - 2*rds*sd*ss + vd + vs) - dl**2*sd*(gm*rdp*sp - rds*ss + sd) + dl*(ad*ss*(gm*rsp*sp + rds*sd - ss) - ad*(gm**2*vp + 2*gm*sp*(rdp*sd - rsp*ss) - 2*rds*sd*ss + vd + vs) + as*sd*(gm*rdp*sp - rds*ss + sd)))) + sD*sS*vS*(omegaD*zetaDP - zPS*zeta**2)*(rDP*sD*sP*(ad*vP*(ad*(gm**2*vp + 2*gm*sp*(rdp*sd - rsp*ss) - 2*rds*sd*ss + vd + vs) + dl*sd*(gm*rdp*sp - rds*ss + sd)) + vD*(gm**2*vp + 2*gm*sp*(rdp*sd - rsp*ss) - 2*rds*sd*ss + vd + vs)) - vD*vP*(2*ad*(gm**2*vp + 2*gm*sp*(rdp*sd - rsp*ss) - 2*rds*sd*ss + vd + vs) + dl*sd*(gm*rdp*sp - rds*ss + sd)))) + sD**2*vP*vD*vP*vS*zeta**2*zetaDD**(13/2)*(as**2*(-gm**2*vp - 2*gm*sp*(rdp*sd - rsp*ss) + 2*rds*sd*ss - vd - vs) + as*dl*(gm**2*vp + 2*gm*sp*(rdp*sd - rsp*ss) - 2*rds*sd*ss + vd + vs - ss*(gm*rsp*sp + rds*sd - ss)) + dl**2*ss*(gm*rsp*sp + rds*sd - ss))*(gD*hS*sqrt(zetaDD)*zetaDS - hS*zSP*zeta*sqrt(zetaDD) + zeta*zetaDD) + sD**2*sP*sS*vD*zetaDD**5*(-gD*zetaDD*(omegaD*rSP*zeta**2*sqrt(zetaDD) + omegaD*zetaDD**(3/2)*zetaSP + zeta**2*zetaDD*(-2*hD*rSP + hP*rDS + hS*rDP)) + zeta**3*zetaDD**(3/2)*(hP*hS + 2*rSP*wS + rSP))*(rSP*sP*sS*(vP*(as**2*(-gm**2*vp - 2*gm*sp*(rdp*sd - rsp*ss) + 2*rds*sd*ss - vd - vs) + as*dl*(gm**2*vp + 2*gm*sp*(rdp*sd - rsp*ss) - 2*rds*sd*ss + vd + vs - ss*(gm*rsp*sp + rds*sd - ss)) + dl**2*ss*(gm*rsp*sp + rds*sd - ss)) - vS*(gm**2*vp + 2*gm*sp*(rdp*sd - rsp*ss) - 2*rds*sd*ss + vd + vs)) + vP*vS*(2*as*(gm**2*vp + 2*gm*sp*(rdp*sd - rsp*ss) - 2*rds*sd*ss + vd + vs) - dl*(gm**2*vp - gm*rsp*sp*ss + 2*gm*sp*(rdp*sd - rsp*ss) - 3*rds*sd*ss + vd + 2*vs))) - sD**2*sS**2*vD*vP*vS*zeta**2*zetaDD**(11/2)*(gD*hP*zetaDD**(3/2)*zetaDP - hS*zSP*zeta*zetaDD**(3/2) + zeta*zetaDD**2 + zeta*zetaDD*(-hP**2 + hS**2))*(gm**2*vp + 2*gm*sp*(rdp*sd - rsp*ss) - 2*rds*sd*ss + vd + vs)))/(dl**3*sD**2*vP*sS**2*vD*vP*vS*zeta**3*zetaDD**(15/2)*zetaSS**(15/2))
  pbeta_s1 <- sweep(xs, MARGIN = 1, (LD*zetaDD**2*zetaSS**(3/2)*(ad*sP*sS*(gS*zetaDS - zDP*zeta) - gS*sD*sP*zetaSS*(as - dl) + sD*sS*(gS*zetaSP - zPD*zeta)) + LS*zetaDD**(3/2)*zetaSS**2*(-ad*gD*sP*sS*zetaDD + sD*sP*(as - dl)*(gD*zetaDS - zSP*zeta) + sD*sS*(gD*zetaDP - zPS*zeta)))/(dl*sD*sP*sS*zeta*zetaDD**2*zetaSS**2), `*`)
  pgamma <- (LD*zetaDD**(15/2)*(ad*dl*sD*vP*sS**2*vD*vP*vS*zeta**2*zetaSS**7*(gS*zetaDS - zDP*zeta)*(LP*gm + xbd + xbp*gm - xbs - dl*(LP + xbp)) - ad*vP*sS**2*vD*vP*vS*zeta**2*zetaSS**(13/2)*(-ad*(gm**2*vp + 2*gm*sp*(rdp*sd - rsp*ss) - 2*rds*sd*ss + vd + vs) + dl**2*rdp*sd*sp + dl*(ad*sp*(gm*sp + rdp*sd - rsp*ss) - sd*(gm*rdp*sp - rds*ss + sd)))*(gS*hD*zetaDS*sqrt(zetaSS) - hD*zDP*zeta*sqrt(zetaSS) + zeta*zetaSS) + dl*sD**2*sP*sS**2*vD*vP*vS*zeta**2*zetaSS**7*(gS*zetaSP - zPD*zeta)*(LP*gm + xbd + xbp*gm - xbs - dl*(LP + xbp)) + gS*sD*sP*zetaSS**(15/2)*(-as*dl*sD*sP*sS*vD*vP*vS*zeta**2*sqrt(zetaSS)*(LP*gm + xbd + xbp*gm - xbs - dl*(LP + xbp)) + as*hS*sD*sP*vD*vP*vS*zeta**2*sqrt(zetaSS)*(-as*(gm**2*vp + 2*gm*sp*(rdp*sd - rsp*ss) - 2*rds*sd*ss + vd + vs) + dl**2*rsp*sp*ss + dl*(as*sp*(gm*sp + rdp*sd - rsp*ss) - ss*(gm*rsp*sp + rds*sd - ss))) + sD*sS*vD*(omegaS*zetaSP - zPD*zeta**2)*(rSP*sP*sS*(as*vP*(-as*(gm**2*vp + 2*gm*sp*(rdp*sd - rsp*ss) - 2*rds*sd*ss + vd + vs) + dl**2*rsp*sp*ss + dl*(as*sp*(gm*sp + rdp*sd - rsp*ss) - ss*(gm*rsp*sp + rds*sd - ss))) - vS*(-dl*sp*(gm*sp + rdp*sd - rsp*ss) + gm**2*vp + 2*gm*sp*(rdp*sd - rsp*ss) - 2*rds*sd*ss + vd + vs)) - vP*vS*(-2*as*(gm**2*vp + 2*gm*sp*(rdp*sd - rsp*ss) - 2*rds*sd*ss + vd + vs) + dl**2*rsp*sp*ss + dl*(2*as*sp*(gm*sp + rdp*sd - rsp*ss) - gm*rsp*sp*ss - rds*sd*ss + vs))) + sP*sS*vP*(omegaS*zetaDS - zDP*zeta**2)*(rDS*sD*sS*(ad*vS*(-ad*(gm**2*vp + 2*gm*sp*(rdp*sd - rsp*ss) - 2*rds*sd*ss + vd + vs) + dl**2*rdp*sd*sp + dl*(ad*sp*(gm*sp + rdp*sd - rsp*ss) - sd*(gm*rdp*sp - rds*ss + sd))) + as*vD*(-as*(gm**2*vp + 2*gm*sp*(rdp*sd - rsp*ss) - 2*rds*sd*ss + vd + vs) + dl**2*rsp*sp*ss + dl*(as*sp*(gm*sp + rdp*sd - rsp*ss) - ss*(gm*rsp*sp + rds*sd - ss)))) + vD*vS*(2*ad*as*(gm**2*vp + 2*gm*sp*(rdp*sd - rsp*ss) - 2*rds*sd*ss + vd + vs) - dl**2*sp*(ad*rsp*ss + as*rdp*sd) + dl*(-2*ad*as*sp*(gm*sp + rdp*sd - rsp*ss) + ad*ss*(gm*rsp*sp + rds*sd - ss) + as*sd*(gm*rdp*sp - rds*ss + sd))))) + sD**2*sS**2*vD*vP*vS*zeta**2*zetaSS**(11/2)*(gS*hP*zetaSP*zetaSS**(3/2) - hD*zDP*zeta*zetaSS**(3/2) + zeta*zetaSS**2 + zeta*zetaSS*(hD**2 - hP**2))*(-dl*sp*(gm*sp + rdp*sd - rsp*ss) + gm**2*vp + 2*gm*sp*(rdp*sd - rsp*ss) - 2*rds*sd*ss + vd + vs) - sD*sP*sS**2*vS*zetaSS**5*(-gS*zetaSS*(omegaS*rDP*zeta**2*sqrt(zetaSS) + omegaS*zetaDP*zetaSS**(3/2) + zeta**2*zetaSS*(hD*rSP + hP*rDS - 2*hS*rDP)) + zeta**3*zetaSS**(3/2)*(hD*hP + 2*rDP*wD + rDP))*(rDP*sD*sP*(ad*vP*(-ad*(gm**2*vp + 2*gm*sp*(rdp*sd - rsp*ss) - 2*rds*sd*ss + vd + vs) + dl**2*rdp*sd*sp + dl*(ad*sp*(gm*sp + rdp*sd - rsp*ss) - sd*(gm*rdp*sp - rds*ss + sd))) - vD*(-dl*sp*(gm*sp + rdp*sd - rsp*ss) + gm**2*vp + 2*gm*sp*(rdp*sd - rsp*ss) - 2*rds*sd*ss + vd + vs)) - vD*vP*(-2*ad*(gm**2*vp + 2*gm*sp*(rdp*sd - rsp*ss) - 2*rds*sd*ss + vd + vs) + dl**2*rdp*sd*sp + dl*(2*ad*sp*(gm*sp + rdp*sd - rsp*ss) - gm*rdp*sd*sp + rds*sd*ss - vd)))) + LS*zetaSS**(15/2)*(as*dl*sD**2*vP*sS*vD*vP*vS*zeta**2*zetaDD**7*(gD*zetaDS - zSP*zeta)*(LP*gm + xbd + xbp*gm - xbs - dl*(LP + xbp)) - as*sD**2*vP*vD*vP*vS*zeta**2*zetaDD**(13/2)*(-as*(gm**2*vp + 2*gm*sp*(rdp*sd - rsp*ss) - 2*rds*sd*ss + vd + vs) + dl**2*rsp*sp*ss + dl*(as*sp*(gm*sp + rdp*sd - rsp*ss) - ss*(gm*rsp*sp + rds*sd - ss)))*(gD*hS*sqrt(zetaDD)*zetaDS - hS*zSP*zeta*sqrt(zetaDD) + zeta*zetaDD) + dl*sD**2*sP*sS**2*vD*vP*vS*zeta**2*zetaDD**7*(gD*zetaDP - zPS*zeta)*(LP*gm + xbd + xbp*gm - xbs - dl*(LP + xbp)) + gD*sP*sS*zetaDD**(15/2)*(-ad*dl*sD*sP*sS*vD*vP*vS*zeta**2*sqrt(zetaDD)*(LP*gm + xbd + xbp*gm - xbs - dl*(LP + xbp)) + ad*hD*sP*sS*vD*vP*vS*zeta**2*sqrt(zetaDD)*(-ad*(gm**2*vp + 2*gm*sp*(rdp*sd - rsp*ss) - 2*rds*sd*ss + vd + vs) + dl**2*rdp*sd*sp + dl*(ad*sp*(gm*sp + rdp*sd - rsp*ss) - sd*(gm*rdp*sp - rds*ss + sd))) + sD*sP*vP*(omegaD*zetaDS - zSP*zeta**2)*(rDS*sD*sS*(ad*vS*(-ad*(gm**2*vp + 2*gm*sp*(rdp*sd - rsp*ss) - 2*rds*sd*ss + vd + vs) + dl**2*rdp*sd*sp + dl*(ad*sp*(gm*sp + rdp*sd - rsp*ss) - sd*(gm*rdp*sp - rds*ss + sd))) + as*vD*(-as*(gm**2*vp + 2*gm*sp*(rdp*sd - rsp*ss) - 2*rds*sd*ss + vd + vs) + dl**2*rsp*sp*ss + dl*(as*sp*(gm*sp + rdp*sd - rsp*ss) - ss*(gm*rsp*sp + rds*sd - ss)))) + vD*vS*(2*ad*as*(gm**2*vp + 2*gm*sp*(rdp*sd - rsp*ss) - 2*rds*sd*ss + vd + vs) - dl**2*sp*(ad*rsp*ss + as*rdp*sd) + dl*(-2*ad*as*sp*(gm*sp + rdp*sd - rsp*ss) + ad*ss*(gm*rsp*sp + rds*sd - ss) + as*sd*(gm*rdp*sp - rds*ss + sd)))) + sD*sS*vS*(omegaD*zetaDP - zPS*zeta**2)*(rDP*sD*sP*(ad*vP*(-ad*(gm**2*vp + 2*gm*sp*(rdp*sd - rsp*ss) - 2*rds*sd*ss + vd + vs) + dl**2*rdp*sd*sp + dl*(ad*sp*(gm*sp + rdp*sd - rsp*ss) - sd*(gm*rdp*sp - rds*ss + sd))) - vD*(-dl*sp*(gm*sp + rdp*sd - rsp*ss) + gm**2*vp + 2*gm*sp*(rdp*sd - rsp*ss) - 2*rds*sd*ss + vd + vs)) - vD*vP*(-2*ad*(gm**2*vp + 2*gm*sp*(rdp*sd - rsp*ss) - 2*rds*sd*ss + vd + vs) + dl**2*rdp*sd*sp + dl*(2*ad*sp*(gm*sp + rdp*sd - rsp*ss) - gm*rdp*sd*sp + rds*sd*ss - vd)))) - sD**2*sP*sS*vD*zetaDD**5*(-gD*zetaDD*(omegaD*rSP*zeta**2*sqrt(zetaDD) + omegaD*zetaDD**(3/2)*zetaSP + zeta**2*zetaDD*(-2*hD*rSP + hP*rDS + hS*rDP)) + zeta**3*zetaDD**(3/2)*(hP*hS + 2*rSP*wS + rSP))*(rSP*sP*sS*(as*vP*(-as*(gm**2*vp + 2*gm*sp*(rdp*sd - rsp*ss) - 2*rds*sd*ss + vd + vs) + dl**2*rsp*sp*ss + dl*(as*sp*(gm*sp + rdp*sd - rsp*ss) - ss*(gm*rsp*sp + rds*sd - ss))) - vS*(-dl*sp*(gm*sp + rdp*sd - rsp*ss) + gm**2*vp + 2*gm*sp*(rdp*sd - rsp*ss) - 2*rds*sd*ss + vd + vs)) - vP*vS*(-2*as*(gm**2*vp + 2*gm*sp*(rdp*sd - rsp*ss) - 2*rds*sd*ss + vd + vs) + dl**2*rsp*sp*ss + dl*(2*as*sp*(gm*sp + rdp*sd - rsp*ss) - gm*rsp*sp*ss - rds*sd*ss + vs))) + sD**2*sS**2*vD*vP*vS*zeta**2*zetaDD**(11/2)*(gD*hP*zetaDD**(3/2)*zetaDP - hS*zSP*zeta*zetaDD**(3/2) + zeta*zetaDD**2 + zeta*zetaDD*(-hP**2 + hS**2))*(-dl*sp*(gm*sp + rdp*sd - rsp*ss) + gm**2*vp + 2*gm*sp*(rdp*sd - rsp*ss) - 2*rds*sd*ss + vd + vs)))/(dl**3*sD**2*vP*sS**2*vD*vP*vS*zeta**3*zetaDD**(15/2)*zetaSS**(15/2))
  pbeta_p1 <- sweep(-xp, MARGIN = 1, gm*(LD*zetaDD**2*zetaSS**(3/2)*(ad*sP*sS*(gS*zetaDS - zDP*zeta) - as*gS*sD*sP*zetaSS + sD*sS*(gS*zetaSP - zPD*zeta)) + LS*zetaDD**(3/2)*zetaSS**2*(-ad*gD*sP*sS*zetaDD + as*sD*sP*(gD*zetaDS - zSP*zeta) + sD*sS*(gD*zetaDP - zPS*zeta)))/(dl*sD*sP*sS*zeta*zetaDD**2*zetaSS**2), `*`)
  pvar_d <- -(-LD*zetaDD**(11/2)*(gS*sD*sP*zetaSS**(11/2)*(as*hS*sD*sP*vD*vP*vS*zeta**2*sqrt(zetaSS)*(as*(gm*rdp*sp - rds*ss + sd) + dl*rds*ss) + sD*sS*vD*(omegaS*zetaSP - zPD*zeta**2)*(rSP*sP*sS*(as*vP*(as*(gm*rdp*sp - rds*ss + sd) + dl*rds*ss) + vS*(gm*rdp*sp - rds*ss + sd)) - vP*vS*(2*as*(gm*rdp*sp - rds*ss + sd) + dl*rds*ss)) + sP*sS*vP*(omegaS*zetaDS - zDP*zeta**2)*(rDS*sD*sS*(as*vD*(as*(gm*rdp*sp - rds*ss + sd) + dl*rds*ss) + vS*(ad**2*(gm*rdp*sp - rds*ss + sd) + ad*dl*(gm*rdp*sp - rds*ss + 2*sd) + dl**2*sd)) - vD*vS*(2*ad*as*(gm*rdp*sp - rds*ss + sd) + dl**2*rds*ss + dl*(ad*rds*ss + as*(gm*rdp*sp - rds*ss + 2*sd))))) - sD**2*sS**2*vD*vP*vS*zeta**2*zetaSS**(7/2)*(gm*rdp*sp - rds*ss + sd)*(gS*hP*zetaSP*zetaSS**(3/2) - hD*zDP*zeta*zetaSS**(3/2) + zeta*zetaSS**2 + zeta*zetaSS*(hD**2 - hP**2)) + sD*sP*sS**2*vS*zetaSS**3*(gS*zetaSS*(omegaS*rDP*zeta**2*sqrt(zetaSS) + omegaS*zetaDP*zetaSS**(3/2) + zeta**2*zetaSS*(hD*rSP + hP*rDS - 2*hS*rDP)) - zeta**3*zetaSS**(3/2)*(hD*hP + 2*rDP*wD + rDP))*(rDP*sD*sP*(vD*(gm*rdp*sp - rds*ss + sd) + vP*(ad**2*(gm*rdp*sp - rds*ss + sd) + ad*dl*(gm*rdp*sp - rds*ss + 2*sd) + dl**2*sd)) - vD*vP*(2*ad*(gm*rdp*sp - rds*ss + sd) + dl*(gm*rdp*sp - rds*ss + 2*sd))) - vP*sS**2*vD*vP*vS*zeta**2*zetaSS**(9/2)*(ad**2*(gm*rdp*sp - rds*ss + sd) + ad*dl*(gm*rdp*sp - rds*ss + 2*sd) + dl**2*sd)*(gS*hD*zetaDS*sqrt(zetaSS) - hD*zDP*zeta*sqrt(zetaSS) + zeta*zetaSS)) + LS*zetaSS**(11/2)*(as*sD**2*vP*vD*vP*vS*zeta**2*zetaDD**(9/2)*(as*(gm*rdp*sp - rds*ss + sd) + dl*rds*ss)*(gD*hS*sqrt(zetaDD)*zetaDS - hS*zSP*zeta*sqrt(zetaDD) + zeta*zetaDD) - gD*sP*sS*zetaDD**(11/2)*(hD*sP*sS*vD*vP*vS*zeta**2*sqrt(zetaDD)*(ad**2*(gm*rdp*sp - rds*ss + sd) + ad*dl*(gm*rdp*sp - rds*ss + 2*sd) + dl**2*sd) + sD*sP*vP*(omegaD*zetaDS - zSP*zeta**2)*(rDS*sD*sS*(as*vD*(as*(gm*rdp*sp - rds*ss + sd) + dl*rds*ss) + vS*(ad**2*(gm*rdp*sp - rds*ss + sd) + ad*dl*(gm*rdp*sp - rds*ss + 2*sd) + dl**2*sd)) - vD*vS*(2*ad*as*(gm*rdp*sp - rds*ss + sd) + dl**2*rds*ss + dl*(ad*rds*ss + as*(gm*rdp*sp - rds*ss + 2*sd)))) + sD*sS*vS*(omegaD*zetaDP - zPS*zeta**2)*(rDP*sD*sP*(vD*(gm*rdp*sp - rds*ss + sd) + vP*(ad**2*(gm*rdp*sp - rds*ss + sd) + ad*dl*(gm*rdp*sp - rds*ss + 2*sd) + dl**2*sd)) - vD*vP*(2*ad*(gm*rdp*sp - rds*ss + sd) + dl*(gm*rdp*sp - rds*ss + 2*sd)))) - sD**2*sP*sS*vD*zetaDD**3*(gD*zetaDD*(omegaD*rSP*zeta**2*sqrt(zetaDD) + omegaD*zetaDD**(3/2)*zetaSP + zeta**2*zetaDD*(-2*hD*rSP + hP*rDS + hS*rDP)) - zeta**3*zetaDD**(3/2)*(hP*hS + 2*rSP*wS + rSP))*(rSP*sP*sS*(as*vP*(as*(gm*rdp*sp - rds*ss + sd) + dl*rds*ss) + vS*(gm*rdp*sp - rds*ss + sd)) - vP*vS*(2*as*(gm*rdp*sp - rds*ss + sd) + dl*rds*ss)) + sD**2*sS**2*vD*vP*vS*zeta**2*zetaDD**(7/2)*(gm*rdp*sp - rds*ss + sd)*(gD*hP*zetaDD**(3/2)*zetaDP - hS*zSP*zeta*zetaDD**(3/2) + zeta*zetaDD**2 - zeta*zetaDD*(hP**2 - hS**2))))/(2*dl**2*sD**2*vP*sS**2*sd*vD*vP*vS*zeta**3*zetaDD**(11/2)*zetaSS**(11/2))
  pvar_s <- (LD*zetaDD**(11/2)*(ad*vP*sS**2*vD*vP*vS*zeta**2*zetaSS**(9/2)*(ad*(gm*rsp*sp + rds*sd - ss) + dl*rds*sd)*(gS*hD*zetaDS*sqrt(zetaSS) - hD*zDP*zeta*sqrt(zetaSS) + zeta*zetaSS) + gS*sD*sP*zetaSS**(11/2)*(hS*sD*sP*vD*vP*vS*zeta**2*sqrt(zetaSS)*(as**2*(-gm*rsp*sp - rds*sd + ss) + as*dl*(gm*rsp*sp + rds*sd - 2*ss) + dl**2*ss) + sD*sS*vD*(omegaS*zetaSP - zPD*zeta**2)*(rSP*sP*sS*(vP*(as**2*(-gm*rsp*sp - rds*sd + ss) + as*dl*(gm*rsp*sp + rds*sd - 2*ss) + dl**2*ss) - vS*(gm*rsp*sp + rds*sd - ss)) + vP*vS*(2*as*(gm*rsp*sp + rds*sd - ss) - dl*(gm*rsp*sp + rds*sd - 2*ss))) - sP*sS*vP*(omegaS*zetaDS - zDP*zeta**2)*(rDS*sD*sS*(ad*vS*(ad*(gm*rsp*sp + rds*sd - ss) + dl*rds*sd) - vD*(as**2*(-gm*rsp*sp - rds*sd + ss) + as*dl*(gm*rsp*sp + rds*sd - 2*ss) + dl**2*ss)) + vD*vS*(-2*ad*as*(gm*rsp*sp + rds*sd - ss) + dl**2*rds*sd + dl*(ad*(gm*rsp*sp + rds*sd - 2*ss) - as*rds*sd)))) + sD**2*sS**2*vD*vP*vS*zeta**2*zetaSS**(7/2)*(gm*rsp*sp + rds*sd - ss)*(gS*hP*zetaSP*zetaSS**(3/2) - hD*zDP*zeta*zetaSS**(3/2) + zeta*zetaSS**2 + zeta*zetaSS*(hD**2 - hP**2)) + sD*sP*sS**2*vS*zetaSS**3*(-gS*zetaSS*(omegaS*rDP*zeta**2*sqrt(zetaSS) + omegaS*zetaDP*zetaSS**(3/2) + zeta**2*zetaSS*(hD*rSP + hP*rDS - 2*hS*rDP)) + zeta**3*zetaSS**(3/2)*(hD*hP + 2*rDP*wD + rDP))*(rDP*sD*sP*(ad*vP*(ad*(gm*rsp*sp + rds*sd - ss) + dl*rds*sd) + vD*(gm*rsp*sp + rds*sd - ss)) - vD*vP*(2*ad*(gm*rsp*sp + rds*sd - ss) + dl*rds*sd))) - LS*zetaSS**(11/2)*(gD*sP*sS*zetaDD**(11/2)*(ad*hD*sP*sS*vD*vP*vS*zeta**2*sqrt(zetaDD)*(ad*(gm*rsp*sp + rds*sd - ss) + dl*rds*sd) + sD*sP*vP*(omegaD*zetaDS - zSP*zeta**2)*(rDS*sD*sS*(ad*vS*(ad*(gm*rsp*sp + rds*sd - ss) + dl*rds*sd) - vD*(as**2*(-gm*rsp*sp - rds*sd + ss) + as*dl*(gm*rsp*sp + rds*sd - 2*ss) + dl**2*ss)) + vD*vS*(-2*ad*as*(gm*rsp*sp + rds*sd - ss) + dl**2*rds*sd + dl*(ad*(gm*rsp*sp + rds*sd - 2*ss) - as*rds*sd))) + sD*sS*vS*(omegaD*zetaDP - zPS*zeta**2)*(rDP*sD*sP*(ad*vP*(ad*(gm*rsp*sp + rds*sd - ss) + dl*rds*sd) + vD*(gm*rsp*sp + rds*sd - ss)) - vD*vP*(2*ad*(gm*rsp*sp + rds*sd - ss) + dl*rds*sd))) + sD**2*vP*vD*vP*vS*zeta**2*zetaDD**(9/2)*(as**2*(-gm*rsp*sp - rds*sd + ss) + as*dl*(gm*rsp*sp + rds*sd - 2*ss) + dl**2*ss)*(gD*hS*sqrt(zetaDD)*zetaDS - hS*zSP*zeta*sqrt(zetaDD) + zeta*zetaDD) + sD**2*sP*sS*vD*zetaDD**3*(-gD*zetaDD*(omegaD*rSP*zeta**2*sqrt(zetaDD) + omegaD*zetaDD**(3/2)*zetaSP + zeta**2*zetaDD*(-2*hD*rSP + hP*rDS + hS*rDP)) + zeta**3*zetaDD**(3/2)*(hP*hS + 2*rSP*wS + rSP))*(rSP*sP*sS*(vP*(as**2*(-gm*rsp*sp - rds*sd + ss) + as*dl*(gm*rsp*sp + rds*sd - 2*ss) + dl**2*ss) - vS*(gm*rsp*sp + rds*sd - ss)) + vP*vS*(2*as*(gm*rsp*sp + rds*sd - ss) - dl*(gm*rsp*sp + rds*sd - 2*ss))) - sD**2*sS**2*vD*vP*vS*zeta**2*zetaDD**(7/2)*(gm*rsp*sp + rds*sd - ss)*(gD*hP*zetaDD**(3/2)*zetaDP - hS*zSP*zeta*zetaDD**(3/2) + zeta*zetaDD**2 + zeta*zetaDD*(-hP**2 + hS**2))))/(2*dl**2*sD**2*vP*sS**2*ss*vD*vP*vS*zeta**3*zetaDD**(11/2)*zetaSS**(11/2))
  pvar_p <- -gm*(LD*zetaDD**(11/2)*(ad*vP*sS**2*vD*vP*vS*zeta**2*zetaSS**(9/2)*(ad*(gm*sp + rdp*sd - rsp*ss) + dl*rdp*sd)*(gS*hD*zetaDS*sqrt(zetaSS) - hD*zDP*zeta*sqrt(zetaSS) + zeta*zetaSS) - gS*sD*sP*zetaSS**(11/2)*(as*hS*sD*sP*vD*vP*vS*zeta**2*sqrt(zetaSS)*(as*(gm*sp + rdp*sd - rsp*ss) + dl*rsp*ss) + sD*sS*vD*(omegaS*zetaSP - zPD*zeta**2)*(rSP*sP*sS*(as*vP*(as*(gm*sp + rdp*sd - rsp*ss) + dl*rsp*ss) + vS*(gm*sp + rdp*sd - rsp*ss)) - vP*vS*(2*as*(gm*sp + rdp*sd - rsp*ss) + dl*rsp*ss)) + sP*sS*vP*(omegaS*zetaDS - zDP*zeta**2)*(rDS*sD*sS*(ad*vS*(ad*(gm*sp + rdp*sd - rsp*ss) + dl*rdp*sd) + as*vD*(as*(gm*sp + rdp*sd - rsp*ss) + dl*rsp*ss)) - vD*vS*(2*ad*as*(gm*sp + rdp*sd - rsp*ss) + dl*(ad*rsp*ss + as*rdp*sd)))) + sD**2*sS**2*vD*vP*vS*zeta**2*zetaSS**(7/2)*(gm*sp + rdp*sd - rsp*ss)*(gS*hP*zetaSP*zetaSS**(3/2) - hD*zDP*zeta*zetaSS**(3/2) + zeta*zetaSS**2 + zeta*zetaSS*(hD**2 - hP**2)) - sD*sP*sS**2*vS*zetaSS**3*(gS*zetaSS*(omegaS*rDP*zeta**2*sqrt(zetaSS) + omegaS*zetaDP*zetaSS**(3/2) + zeta**2*zetaSS*(hD*rSP + hP*rDS - 2*hS*rDP)) - zeta**3*zetaSS**(3/2)*(hD*hP + 2*rDP*wD + rDP))*(rDP*sD*sP*(ad*vP*(ad*(gm*sp + rdp*sd - rsp*ss) + dl*rdp*sd) + vD*(gm*sp + rdp*sd - rsp*ss)) - vD*vP*(2*ad*(gm*sp + rdp*sd - rsp*ss) + dl*rdp*sd))) + LS*zetaSS**(11/2)*(as*sD**2*vP*vD*vP*vS*zeta**2*zetaDD**(9/2)*(as*(gm*sp + rdp*sd - rsp*ss) + dl*rsp*ss)*(gD*hS*sqrt(zetaDD)*zetaDS - hS*zSP*zeta*sqrt(zetaDD) + zeta*zetaDD) - gD*sP*sS*zetaDD**(11/2)*(ad*hD*sP*sS*vD*vP*vS*zeta**2*sqrt(zetaDD)*(ad*(gm*sp + rdp*sd - rsp*ss) + dl*rdp*sd) + sD*sP*vP*(omegaD*zetaDS - zSP*zeta**2)*(rDS*sD*sS*(ad*vS*(ad*(gm*sp + rdp*sd - rsp*ss) + dl*rdp*sd) + as*vD*(as*(gm*sp + rdp*sd - rsp*ss) + dl*rsp*ss)) - vD*vS*(2*ad*as*(gm*sp + rdp*sd - rsp*ss) + dl*(ad*rsp*ss + as*rdp*sd))) + sD*sS*vS*(omegaD*zetaDP - zPS*zeta**2)*(rDP*sD*sP*(ad*vP*(ad*(gm*sp + rdp*sd - rsp*ss) + dl*rdp*sd) + vD*(gm*sp + rdp*sd - rsp*ss)) - vD*vP*(2*ad*(gm*sp + rdp*sd - rsp*ss) + dl*rdp*sd))) - sD**2*sP*sS*vD*zetaDD**3*(gD*zetaDD*(omegaD*rSP*zeta**2*sqrt(zetaDD) + omegaD*zetaDD**(3/2)*zetaSP + zeta**2*zetaDD*(-2*hD*rSP + hP*rDS + hS*rDP)) - zeta**3*zetaDD**(3/2)*(hP*hS + 2*rSP*wS + rSP))*(rSP*sP*sS*(as*vP*(as*(gm*sp + rdp*sd - rsp*ss) + dl*rsp*ss) + vS*(gm*sp + rdp*sd - rsp*ss)) - vP*vS*(2*as*(gm*sp + rdp*sd - rsp*ss) + dl*rsp*ss)) + sD**2*sS**2*vD*vP*vS*zeta**2*zetaDD**(7/2)*(gm*sp + rdp*sd - rsp*ss)*(gD*hP*zetaDD**(3/2)*zetaDP - hS*zSP*zeta*zetaDD**(3/2) + zeta*zetaDD**2 - zeta*zetaDD*(hP**2 - hS**2))))/(2*dl**2*sD**2*vP*sS**2*sp*vD*vP*vS*zeta**3*zetaDD**(11/2)*zetaSS**(11/2))
  Dl <- cbind(palpha_d, pbeta_d1, palpha_s, pbeta_s1, pgamma, pbeta_p1, pvar_d, pvar_s, pvar_p)
  if (object@correlated_shocks) {
      prho_ds <- sd*ss*(LD*zetaDD**(11/2)*(ad*vP*sS**2*vD*vP*vS*zeta**2*zetaSS**(9/2)*(ad + dl)*(gS*hD*zetaDS*sqrt(zetaSS) - hD*zDP*zeta*sqrt(zetaSS) + zeta*zetaSS) - gS*sD*sP*zetaSS**(11/2)*(as*hS*sD*sP*vD*vP*vS*zeta**2*sqrt(zetaSS)*(as - dl) + sD*sS*vD*(omegaS*zetaSP - zPD*zeta**2)*(rSP*sP*sS*(as*vP*(as - dl) + vS) - vP*vS*(2*as - dl)) + sP*sS*vP*(omegaS*zetaDS - zDP*zeta**2)*(rDS*sD*sS*(ad*vS*(ad + dl) + as*vD*(as - dl)) + vD*vS*(-2*ad*as + dl**2 + dl*(ad - as)))) + sD**2*sS**2*vD*vP*vS*zeta**2*zetaSS**(7/2)*(gS*hP*zetaSP*zetaSS**(3/2) - hD*zDP*zeta*zetaSS**(3/2) + zeta*zetaSS**2 + zeta*zetaSS*(hD**2 - hP**2)) + sD*sP*sS**2*vS*zetaSS**3*(-gS*zetaSS*(omegaS*rDP*zeta**2*sqrt(zetaSS) + omegaS*zetaDP*zetaSS**(3/2) + zeta**2*zetaSS*(hD*rSP + hP*rDS - 2*hS*rDP)) + zeta**3*zetaSS**(3/2)*(hD*hP + 2*rDP*wD + rDP))*(rDP*sD*sP*(ad*vP*(ad + dl) + vD) - vD*vP*(2*ad + dl))) + LS*zetaSS**(11/2)*(as*sD**2*vP*vD*vP*vS*zeta**2*zetaDD**(9/2)*(as - dl)*(gD*hS*sqrt(zetaDD)*zetaDS - hS*zSP*zeta*sqrt(zetaDD) + zeta*zetaDD) - gD*sP*sS*zetaDD**(11/2)*(ad*hD*sP*sS*vD*vP*vS*zeta**2*sqrt(zetaDD)*(ad + dl) + sD*sP*vP*(omegaD*zetaDS - zSP*zeta**2)*(rDS*sD*sS*(ad*vS*(ad + dl) + as*vD*(as - dl)) + vD*vS*(-2*ad*as + dl**2 + dl*(ad - as))) + sD*sS*vS*(omegaD*zetaDP - zPS*zeta**2)*(rDP*sD*sP*(ad*vP*(ad + dl) + vD) - vD*vP*(2*ad + dl))) + sD**2*sP*sS*vD*zetaDD**3*(-gD*zetaDD*(omegaD*rSP*zeta**2*sqrt(zetaDD) + omegaD*zetaDD**(3/2)*zetaSP + zeta**2*zetaDD*(-2*hD*rSP + hP*rDS + hS*rDP)) + zeta**3*zetaDD**(3/2)*(hP*hS + 2*rSP*wS + rSP))*(rSP*sP*sS*(as*vP*(as - dl) + vS) - vP*vS*(2*as - dl)) + sD**2*sS**2*vD*vP*vS*zeta**2*zetaDD**(7/2)*(gD*hP*zetaDD**(3/2)*zetaDP - hS*zSP*zeta*zetaDD**(3/2) + zeta*zetaDD**2 + zeta*zetaDD*(-hP**2 + hS**2))))/(dl**2*sD**2*vP*sS**2*vD*vP*vS*zeta**3*zetaDD**(11/2)*zetaSS**(11/2))
  prho_dp <- -gm*sd*sp*(LD*zetaDD**(11/2)*(ad*vP*sS**2*vD*vP*vS*zeta**2*zetaSS**(9/2)*(ad + dl)*(gS*hD*zetaDS*sqrt(zetaSS) - hD*zDP*zeta*sqrt(zetaSS) + zeta*zetaSS) + gS*sD*sP*zetaSS**(11/2)*(-as**2*hS*sD*sP*vD*vP*vS*zeta**2*sqrt(zetaSS) + sD*sS*vD*(omegaS*zetaSP - zPD*zeta**2)*(2*as*vP*vS - rSP*sP*sS*(as**2*vP + vS)) + sP*sS*vP*(omegaS*zetaDS - zDP*zeta**2)*(as*vD*vS*(2*ad + dl) - rDS*sD*sS*(ad*vS*(ad + dl) + as**2*vD))) + sD**2*sS**2*vD*vP*vS*zeta**2*zetaSS**(7/2)*(gS*hP*zetaSP*zetaSS**(3/2) - hD*zDP*zeta*zetaSS**(3/2) + zeta*zetaSS**2 + zeta*zetaSS*(hD**2 - hP**2)) - sD*sP*sS**2*vS*zetaSS**3*(gS*zetaSS*(omegaS*rDP*zeta**2*sqrt(zetaSS) + omegaS*zetaDP*zetaSS**(3/2) + zeta**2*zetaSS*(hD*rSP + hP*rDS - 2*hS*rDP)) - zeta**3*zetaSS**(3/2)*(hD*hP + 2*rDP*wD + rDP))*(rDP*sD*sP*(ad*vP*(ad + dl) + vD) - vD*vP*(2*ad + dl))) + LS*zetaSS**(11/2)*(as**2*sD**2*vP*vD*vP*vS*zeta**2*zetaDD**(9/2)*(gD*hS*sqrt(zetaDD)*zetaDS - hS*zSP*zeta*sqrt(zetaDD) + zeta*zetaDD) - gD*sP*sS*zetaDD**(11/2)*(ad*hD*sP*sS*vD*vP*vS*zeta**2*sqrt(zetaDD)*(ad + dl) - sD*sP*vP*(omegaD*zetaDS - zSP*zeta**2)*(as*vD*vS*(2*ad + dl) - rDS*sD*sS*(ad*vS*(ad + dl) + as**2*vD)) + sD*sS*vS*(omegaD*zetaDP - zPS*zeta**2)*(rDP*sD*sP*(ad*vP*(ad + dl) + vD) - vD*vP*(2*ad + dl))) + sD**2*sP*sS*vD*zetaDD**3*(2*as*vP*vS - rSP*sP*sS*(as**2*vP + vS))*(gD*zetaDD*(omegaD*rSP*zeta**2*sqrt(zetaDD) + omegaD*zetaDD**(3/2)*zetaSP + zeta**2*zetaDD*(-2*hD*rSP + hP*rDS + hS*rDP)) - zeta**3*zetaDD**(3/2)*(hP*hS + 2*rSP*wS + rSP)) + sD**2*sS**2*vD*vP*vS*zeta**2*zetaDD**(7/2)*(gD*hP*zetaDD**(3/2)*zetaDP - hS*zSP*zeta*zetaDD**(3/2) + zeta*zetaDD**2 - zeta*zetaDD*(hP**2 - hS**2))))/(dl**2*sD**2*vP*sS**2*vD*vP*vS*zeta**3*zetaDD**(11/2)*zetaSS**(11/2))
  prho_sp <- gm*sp*ss*(LD*zetaDD**(11/2)*(ad**2*vP*sS**2*vD*vP*vS*zeta**2*zetaSS**(9/2)*(gS*hD*zetaDS*sqrt(zetaSS) - hD*zDP*zeta*sqrt(zetaSS) + zeta*zetaSS) + gS*sD*sP*zetaSS**(11/2)*(-as*hS*sD*sP*vD*vP*vS*zeta**2*sqrt(zetaSS)*(as - dl) - sD*sS*vD*(omegaS*zetaSP - zPD*zeta**2)*(rSP*sP*sS*(as*vP*(as - dl) + vS) - vP*vS*(2*as - dl)) + sP*sS*vP*(omegaS*zetaDS - zDP*zeta**2)*(ad*vD*vS*(2*as - dl) - rDS*sD*sS*(ad**2*vS + as*vD*(as - dl)))) + sD**2*sS**2*vD*vP*vS*zeta**2*zetaSS**(7/2)*(gS*hP*zetaSP*zetaSS**(3/2) - hD*zDP*zeta*zetaSS**(3/2) + zeta*zetaSS**2 + zeta*zetaSS*(hD**2 - hP**2)) - sD*sP*sS**2*vS*zetaSS**3*(2*ad*vD*vP - rDP*sD*sP*(ad**2*vP + vD))*(-gS*zetaSS*(omegaS*rDP*zeta**2*sqrt(zetaSS) + omegaS*zetaDP*zetaSS**(3/2) + zeta**2*zetaSS*(hD*rSP + hP*rDS - 2*hS*rDP)) + zeta**3*zetaSS**(3/2)*(hD*hP + 2*rDP*wD + rDP))) + LS*zetaSS**(11/2)*(as*sD**2*vP*vD*vP*vS*zeta**2*zetaDD**(9/2)*(as - dl)*(gD*hS*sqrt(zetaDD)*zetaDS - hS*zSP*zeta*sqrt(zetaDD) + zeta*zetaDD) + gD*sP*sS*zetaDD**(11/2)*(-ad**2*hD*sP*sS*vD*vP*vS*zeta**2*sqrt(zetaDD) + sD*sP*vP*(omegaD*zetaDS - zSP*zeta**2)*(ad*vD*vS*(2*as - dl) - rDS*sD*sS*(ad**2*vS + as*vD*(as - dl))) + sD*sS*vS*(omegaD*zetaDP - zPS*zeta**2)*(2*ad*vD*vP - rDP*sD*sP*(ad**2*vP + vD))) + sD**2*sP*sS*vD*zetaDD**3*(-gD*zetaDD*(omegaD*rSP*zeta**2*sqrt(zetaDD) + omegaD*zetaDD**(3/2)*zetaSP + zeta**2*zetaDD*(-2*hD*rSP + hP*rDS + hS*rDP)) + zeta**3*zetaDD**(3/2)*(hP*hS + 2*rSP*wS + rSP))*(rSP*sP*sS*(as*vP*(as - dl) + vS) - vP*vS*(2*as - dl)) + sD**2*sS**2*vD*vP*vS*zeta**2*zetaDD**(7/2)*(gD*hP*zetaDD**(3/2)*zetaDP - hS*zSP*zeta*zetaDD**(3/2) + zeta*zetaDD**2 + zeta*zetaDD*(-hP**2 + hS**2))))/(dl**2*sD**2*vP*sS**2*vD*vP*vS*zeta**3*zetaDD**(11/2)*zetaSS**(11/2))
    Dl <- cbind(Dl, prho_ds, prho_dp, prho_sp)
  }
  # nolint end

  colnames(Dl) <- likelihood_variables(object)

  Dl / c(object@L_D + object@L_S)
})