Attribute VB_Name = "Choleskey_decomposition"
Public Sub randomdata()             '(x, nv, nc)

Dim rdata() As Double

nv = 5             '変数の数
nc = 300           'ケースの数
ReDim x(nv, nv)    'ターゲットの相関行列

For i = 1 To nv
x(i, i) = 1
Next i
x(2, 1) = 0.3
x(3, 1) = 0.4
x(3, 2) = 0.3
x(4, 1) = 0.2
x(4, 2) = 0.55
x(4, 3) = 0.32
x(5, 1) = 0.2
x(5, 2) = 0.45
x(5, 3) = 0.33
x(5, 4) = 0.43


Dim y() As Double

ReDim rdata(nv, nc)
Call makedata(nv, nc, rdata)
N = nv

'  ReDim x(N, N)
  ReDim y(N, N)
' コレスキー分解呼び出し
Call chol_main(x, y, nv)
Dim z() As Double
ReDim z(nv, nc)

'compute d=d*r1.
'任意の相関行列の乱数データの作成
For j = 1 To nc
For i = 1 To nv
For k = 1 To nv
  z(i, j) = z(i, j) + rdata(k, j) * y(k, i)
Next k
Next i
Next j
For i = 1 To nv
For j = 1 To nc
   Cells(j + 4, i) = z(i, j)
Next j
Next i
End Sub
Private Sub chol_main(x, y, N)
  Dim IFAULT As Integer, NDIM As Integer, NULLTY As Integer
  Dim a() As Double, b() As Double
  NDIM = N * (N + 1) / 2
  ReDim a(NDIM), U(NDIM), b(N)
  ' a()に相関行列 下三角行列を１行で入れる。
  iI = 0
  For i = 1 To N
  For j = 1 To i
  iI = iI + 1
  a(iI) = x(i, j)
  Next j
  Next i
  For i = 1 To N
    b(i) = i
  Next i
  Call SUBCHL(a, b, N, U, NULLTY, IFAULT, NDIM, DET)
  If IFAULT > 0 Then
     MsgBox ("error in subchl")
     End
  End If
 ' 結果の　U()をy(,)の上三角行列に入れる。
  iI = 0
  For i = 1 To N
  For j = 1 To i
  iI = iI + 1
  y(j, i) = U(iI)
'  MsgBox U(ii)
  Next j
  Next i
End Sub

Private Sub SUBCHL(a, b, N, U, NULLTY, IFAULT, NDIM, DET)
'C     Fortranプログラムをvba に変換した。
'     http://lib.stat.cmu.edu/modules.php?op=modload&name=PostWrap&file=index&page=apstat/
'     http://lib.stat.cmu.edu/apstat/6
'C     REMARK ASR 44  APPL. STATIST. (1982) VOL. 31, NO. 3
'C
'C     A revised and enhanced version of
'C     ALGORITHM AS 6  APPL. STATIST. (1968) VOL. 17, NO. 2
'C
'C     Given a symmetric matrix of order N as lower triangle in A(),
'C     calculates an upper triangle, U(), such that U'U = the sub-matrix
'C     of A whose rows and columns are specified in the integer array
'C     B().
'C     U() may coincide with A().   A() must be +ve semi-definite.
'C     ETA is set to multiplying factor determining effective zero for
'C     a pivot.
'C     NULLTY is returned as number of effective zero pivots.
'C     IFAULT is returned as 1 if N <= 0, 2 if A() is not +ve semi-
'C     definite, otherwise 0 is returned.
'C
'      DOUBLE PRECISION A(NDIM),U(NDIM),DET
'      INTEGER B(N)
'C
'C     Local variables
'C
      Dim ETA As Double, ONE As Double, ZERO As Double, W As Double, ETA2 As Double
'C
'C     The value of ETA below will depend upon the word length of the
'C     computer being used.
'C
      ETA = 0.00000000000001
      ONE = 1#
      ZERO = 0#
'C
      IFAULT = 1
      If (N <= 0) Then GoTo 90
      IFAULT = 2
      NULLTY = 0
      DET = ONE
      j = 1
      k = 0
      ETA2 = ETA * ETA
      For icol = 1 To N
        IJ = b(icol) * (b(icol) - 1) / 2
        iI = IJ + b(icol)
        x = ETA2 * a(iI)
        L = 0
        For irow = 1 To icol
          KK = b(irow) * (b(irow) + 1) / 2
          k = k + 1
          jj = IJ + b(irow)
          W = a(jj)
          m = j
          For i = 1 To irow
            L = L + 1
            If (i = irow) Then GoTo 20
            W = W - U(L) * U(m)
            m = m + 1
10        Next i
20        If (irow = icol) Then GoTo 50
          If (U(L) = ZERO) Then GoTo 30
          U(k) = W / U(L)
          GoTo 40
30        If (W * W > Abs(x * a(KK))) Then GoTo 90
          U(k) = ZERO
40      Next irow
50      If (Abs(W) <= Abs(ETA * a(KK))) Then GoTo 60
        If (W < ZERO) Then GoTo 90
        U(k) = Sqr(W)
        GoTo 70
60      U(k) = ZERO
        NULLTY = NULLTY + 1
70      j = j + icol
        DET = DET * U(k) * U(k)
80    Next icol
'C
      IFAULT = 0
90    Exit Sub
      End Sub

Sub makedata(nv, nc, x)
 '   nv = N
 '   nc = nc
    Randomize
    For j = 1 To nv
      For i = 1 To nc Step 2
      Call rnd_generator(x1, x2)  '******
      x(j, i) = x1
      If i + 1 > nc Then Exit For
      x(j, i + 1) = x2
      Next i
    Next j
End Sub
    
Private Sub rnd_generator(x1, x2)
'極座標法の Marsaglia の改良 Box-Muller よりも速い
'伏見正則『乱数』東京大学出版会, 1989, p65-67
Dim u1 As Double, u2 As Double, v1, v2, s
next1:
u1 = Rnd(1)
u2 = Rnd(1)
v1 = 2 * u1 - 1
v2 = 2 * u2 - 1
s = v1 ^ 2 + v2 ^ 2
If s >= 1 Then GoTo next1
x1 = v1 * Sqr(-2 * Log(s) / s)
x2 = v2 * Sqr(-2 * Log(s) / s)
End Sub