Q: I have an estimate of a variance-covariance matrix and I need the Cholesky decomposition of it. Does anybody have any suggestions on how to program this in GAMS?

There are multiple ways to calculate the Cholesky factorization in GAMS: implement the algorithm in GAMS, use a model and solver, or use the cholesky tool. These three methods are demonstrated in the model below. Thanks to Oyvind Hoveid and Arne Drud for contributing to this model.

Sets
     j          'variables'                /j1*j10/
     n          'observations'             /n1*n15/
     b_l(j,j)   'lower triangular marix';

Alias (j,jj,jjj);

b_l(j,jj) $ (ord(j) >= ord(jj)) = yes;

Parameters
 jn1(j,n)   'variable observations'
 b1(j,jj)   'covariance matrix';

* forming a covariance matrix
jn1(j,n) = normal(0,1);
b1(j,j)  = sum(n, jn1(j,n))/card(n);
b1(j,jj) = sum(n, (jn1(j,n) - b1(j,j))*(jn1(jj,n) - b1(jj,jj)))/card(n);

* Symbols for the Cholesky decomposition algorithm
Sets
    l(j,j)  'defines the lower triangle'
    r(j)    'the pivot row'                  / #j /
    npr(j)  'rows not yet used for pivot'    / #j /
    npc(j)  'columns not yet used for pivot' / #j /
parameters
    w(j,j)  'factor of b'
    ww(j,j) 'part of b not yet factored'
    big     'the diagonal pivot at each point in the factorization'
    piv     'the square root of the pivot'
    tol     'pivot tolerance' / 1e-10 /;

* Cholesky decomposition algorithm in GAMS       
ww(j,jj) = b1(j,jj);
loop(r,
   big = ww(r,r);
   if (big >= tol,
      piv = sqrt(big);
      l(npr,r) = yes;
      npr(r)   = no;
      npc(r)   = no;
      w(r  ,r) = piv;
      w(npr,r) = ww(npr,r) / piv;
      ww(npr,npc) = ww(npr,npc) - ww(npr,r)*ww(npc,r)/big;
      ww(npr,r) = 0;
      ww(r,npc) = 0;
      ww(r,r)   = 0;
   );
);
display w;

* Using a CNS model to determine the Cholesky factors
Variables
   BB(j,jj) 'Cholesky factors'
Equations
   b1_eq(j,jj) 'Cholesky factorization';
 
b1_eq(b_l(j,jj))..
  b1(j,jj) =e= sum(jjj, (BB(j,jjj)$b_l(j,jjj))*(BB(jj,jjj)$b_l(jj,jjj)));

Model cholesky / b1_eq / ;
BB.l(j,j) = 1;

option limrow=0, limcol=0, solprint=silent;
Solve cholesky using cns;

display bb.l;

* Using a cholesky tool to determine the Cholesky factors
execute_unload 'a.gdx', j, b1;
execute 'cholesky a.gdx j b1 b.gdx w';
execute_load 'b.gdx', w;
display w;