Conditional random sampling: annual to monthly

I have a bunch of annual numbers for X[, 1], X[, 2], ..., X[, numberAssetClasses].

I want to generate random identically distributed conditional multi-normal monthly numbers Y that (summed) match the annual numbers (of each year) but are based on an unconditional distribution.

I think my trouble stems from not knowing how to set the unconditional covariance between the monthly Y[, i]s, derived from the covariance matrix for the annual Xs. (Marked "This seems suspect".)

Here is my code. I hope it would yield similar correlation matrices. It doesn't. Any ideas?

library(mvtnorm)

numberPeriods <- 1000
numberAssetClasses <- 5
numberSubs <- 12

set.seed(42)

X <- matrix(rnorm(numberAssetClasses*numberPeriods, 10, 15*sqrt(12)),
           nrow = numberPeriods)

meanX <- c(rep(0, numberAssetClasses))
for (i in 1:numberAssetClasses){
  meanX[i] <- mean(X[, i])
}

covX <- cov(X)

covX

meanY <- c(rep(0, (numberSubs-1)*numberAssetClasses))
for (i in 1:numberAssetClasses){
  meanY[((i-1)*(numberSubs-1)+1):((i)*(numberSubs-1))] <-
    meanX[i]/numberSubs
}

covY <- covX
for (i in 1:numberAssetClasses){
  for (j in 1:numberAssetClasses){
    if (i == j){
      covY[i, j] <- covX[i, j] / numberSubs
    } else {
      covY[i, j] <- covX[i, j] / numberSubs # This seems suspect
    }
  }
}

bigSigma <- matrix(0.0, nrow = (numberAssetClasses)*(numberSubs-1),
                   ncol = (numberAssetClasses)*(numberSubs-1))
for (i in 1:numberAssetClasses){
  for (j in 1:numberAssetClasses){
    if (i == j){
      bigSigma[((i-1)*(numberSubs-1)+1):((i)*(numberSubs-1)),
               ((j-1)*(numberSubs-1)+1):((j)*(numberSubs-1))] <-
        covY[i, j]*(diag(numberSubs-1)-matrix(1/numberSubs, numberSubs-1, numberSubs-1))
    } else {
      bigSigma[((i-1)*(numberSubs-1)+1):((i)*(numberSubs-1)),
               ((j-1)*(numberSubs-1)+1):((j)*(numberSubs-1))] <-
        covY[i, j]*diag(numberSubs-1)      
    }
  }
}

Y <- matrix(0.0, nrow = numberSubs*numberPeriods, ncol = numberAssetClasses)
for (m in 1:numberPeriods){
  for (n in 1:(numberSubs-1)){
    Y[((m-1)*numberSubs+1):((m)*numberSubs-1),
      (1:numberAssetClasses)] <-
      matrix(rmvnorm(1, meanY, bigSigma), ncol = numberAssetClasses)
  }
  for (i in 1:numberAssetClasses){
    Y[m*numberSubs, i] <- numberSubs*meanY[(i-1)*(numberSubs-1)+1] -
      sum(Y[((m-1)*numberSubs+1):((m)*numberSubs-1), i])
  }
}

cor(X)
cor(Y)