ALS - GDP

Hi to every one, nice to know that forums such as this exist. I would need some help with the code witch was written by one professor who is not available. The problem is that I am completely new to GAUSS and don't know how  to input my data properly to run the code. So far I spent a lot of time and effort in vain. The data is a time series of GDP in excel table. I would really appreciated help about this. Thank you. The code is:

proc(12) = als(x, y, tol, progflag);

@ Computes Adaptive Least Squares filter and smoother
for model y = x*beta(t) + eps @
@ for details, see J. Huston McCulloch, "The Kalman Foundations of
Adaptive Least Squares," Sept. 2004, at
 @

@ Finds signal/noise variance ratio rho by ML @
@ x is nXk regressor matrix.  First column = ones for constant term @
@ y is nX1 dependent variable vector @
@ tol is tolerance for precision of log likelihood.  Try .001 or less @
@ progflag = 1 to show ML progress onscreen,  = 0 for nothing onscreen @
@ returns {bf, pf, bs, ps, u, s2eps, rho, serho, tinf, setinf, logl, lr} @
@ bf is nXk matrix whose t-th row is b(t)', the
filter estimate of beta(t) conditional on y(1) ... y(t),
for t = k, ... n.  First k-1 rows are "missing" code.  @
@ pf is nXk^2 super-matrix of coefficient variance matrices.
reshape(pf[t,.],k,k) is covariance of b(t), for
t = k, ... n.  First k-1 rows are "missing" code.  @
@ bs, ps are similarly formatted smoother estimate and covariances. @
@ u is nx1 vector of variance-equalized forecasting errors.
Should be NID(0,s2eps) if model correctly specified. @
@ s2eps is estimated variance of eps @
@ rho is estimated signal/noise variance ratio. @
@ serho = se(rho) @
@ tinf is asymptotic effective sample size.  Asymp. gain = 1/tinf. @
@ setinf = se(tinf) @
@ logl is maximized log likelihood @
@ lr is 2*(logl(estimated rho) - logl(rho = 0))
This is informative, but distribution is not chi-square. @

   local k, n, rhohi, bf, pf, bs, ps, s2eps, logl0, logl, lr, rho, u, tinf ;
   local serho, setinf, d2logl;
   declare matrix alxglob;
   declare matrix alyglob;

   k = cols(x);
   n = rows(x);

   if n /= rows(y);
      format /rdn 10,0;
      print "rows(x) =/ rows(y) in ALS";
      print "rows(x) " rows(x) " rows(y) " rows(y);
      stop;
   endif;

   if n < k+1;
      format /rdn 10,0;
      print "rows(x) < cols(x) + 1 in ALS";
      print "rows(x) " rows(x) " cols(x) " cols(x);
      stop;
   endif;

   rhohi = 1/rows(y);  @ makes T about sqrt(n) @ 
   alxglob = x;
   alyglob = y;

   {rho, logl, d2logl} = hillright(&alsloglik, 0, rhohi, tol, progflag);

   {bf, pf, bs, ps, u, s2eps, logl0} = alsrho(x,y,0,0);   @ that's a zero in logl0 @

   {bf, pf, bs, ps, u, s2eps, logl}  = alsrho(x,y,rho,1);

   lr = 2*(logl - logl0);
   tinf = .5 + sqrt(.25 + 1/rho);
   serho = sqrt(-1/d2logl);
   setinf = serho/((1+4/rho)^.5 * rho^2);

retp(bf, pf, bs, ps, u, s2eps, rho, serho, tinf, setinf, logl, lr);
endp;

 

proc(1) = alsloglik (rho);
@ Computes logl from rho for proc als to use with hillright.
Receives x and y from als as globals alxglob, alyglob.  @

   local bf, pf, bs, ps, s2eps, logl, u;

   {bf, pf, bs, ps, u, s2eps, logl} = alsrho(alxglob, alyglob, rho, 0);

retp(logl);
endp;

 

proc(7) = alsrho(x, y, rho, smooflag);

@ Computes Adaptive Learning filter and/or smoother given rho,
for model y = x*beta(t) + eps @
@ x is nXk regressor matrix.  First column = ones for constant term @
@ y is nX1 dependent variable vector @
@ rho is signal/noise variance ratio @
@ smooflag = 1 to compute smoother,  = 0 else. @
@ returns {bf, pf, bs, ps, u, s2eps, logl} @
@ bf is nXk matrix whose t-th row is b(t), the
filter estimate of beta(t) conditional on y(1) ... y(t),
for t = k, ... n.  First k-1 rows are "missing" code.  @
@ pf is nXk^2 super-matrix of coefficient variance matices.
reshape(pf[t,.],k,k) is covariance of b(t), for
t = k, ... n.  First k-1 rows are "missing" code.  @
@ bs, ps are similarly formatted smoother estimate and covariances,
if smooflag = 1.  Null vectors returned if smooflag = 0. @
@ u is nx1 vector of variance-equalized forecast errors.
First k values are "missing" code. Should be NID(0, s2eps)
if model correct.  @
@ s2eps is estimated variance of eps @
@ logl is log likelihood @

   local k, n, tt, z, w, logl, s2eps, bf, pf, bs, ps, fade, s, wi, missing, u;
   local zt, wt, zstar, wstar, ws, wsi, wiwrtw, tta, t;

   missing = {.};
   k = cols(x);
   n = rows(x);

   if n /= rows(y);
      format /rdn 10,0;
      print "rows(x) =/ rows(y)";
      print "rows(x) " rows(x) " rows(y) " rows(y);
      stop;
   endif;

   if n < k+1;
      format /rdn 10,0;
      print "rows(x) < cols(x) + 1";
      print "rows(x) " rows(x) " cols(x) " cols(x);
      stop;
   endif;

   z = zeros(k,1);
   zt = {};

   w = zeros(k,k);
   wt = {};
   logl = 0;
   s2eps = 0;
   bf = missing * zeros(n,k);
   pf = missing * zeros(n,k*k);
   u = missing * zeros(n,1);
   bs = {};
   ps = {};
   tt = 0;
   tta = {};

   for t(1,n,1);
      fade = (1+rho*tt);
      if t > k;
         s = sqrt(fade*x[t,.]*wi*x[t,.]' + 1);
         logl = logl - ln(s);
         u[t] = (y[t] - x[t,.]*bf[t-1,.]')/s;
      endif;

      tt = tt/fade + 1;
      tta = tta|tt;
      z  =  z/fade + x[t,.]'*y[t];
      zt = zt~z;

      w  =  w/fade + x[t,.]'*x[t,.];
      wt = wt~vecr(w);
      if t >= k;
         wi = inv(w);
         bf[t,.] = (wi*z)';
         pf[t,.]=vecr(wi)';
      endif;
   endfor;

   s2eps = (u[k+1:n]' * u[k+1:n]) / (n-k);
   logl = logl - ((n-k)/2) * (1 + ln(2*pi*s2eps));
   pf = pf * s2eps;

   if smooflag;
      bs = missing*zeros(n,k);
      ps = missing*zeros(n,k*k);
      zstar = zeros(k,1);
      wstar = zeros(k,k);
      t = n;
      do until t < k;
         w = reshape(wt[.,t],k,k);
         wiwrtw = w*inv(w + rho*tta[t]*wstar);
         ws = w + wiwrtw*wstar;
         wsi = inv(ws);
         ps[t,.] = s2eps*vecr(wsi)';
         bs[t,.] = (wsi*(zt[.,t] + wiwrtw*zstar))';
         zstar = wiwrtw*zstar + x[t,.]'*y[t];
         wstar = wiwrtw*wstar + x[t,.]'*x[t,.];
         t = t - 1;
      endo;
   endif;

retp(bf, pf, bs, ps, u, s2eps, logl);
endp;

 

proc(3) = hillright(&f, xlo, xhi, tol, progflag);

@ maximizes f(x) by golden ratio search @
@ loosely after Numerical Recipes. @
@ returns xmax, f(xmax), d2y
d2y is 2nd derivative at xmax
should be negative if xmax > xlo
= missing code if xmax = xlo @
@ may expand to right above xhi but not to left below xlo @
@ runs until dy < tol (Different than hill, hillur!)  @
@ progflag = 1 for onscreen progress, 0 else.  @

   local f:proc, a, x1, x2, x3, x4, dx, y1, y2, y3, y4, dy;
   local dx2, dx3, dy2, dy3, d2y, missing;

   a = (3 - sqrt(5))/2;         @ a = .382... = 1 - golden ratio  @
   x1 = xlo;
   x4 = xhi;
   dx = x4 - x1;
   x2 = x1 + a * dx;
   x3 = x4 - a * dx;
   y1 = f(x1); y2 = f(x2); y3 = f(x3); y4 = f(x4);
   dy = tol + 1000;

   if progflag;
      format /rdn 13,6;
      print;
      "Maximization progress";
   endif;

   right:

   if y1 < y4 and y2 < y4 and y3 < y4;
      x2 = x3; y2 = y3;
      x3 = x4; y3 = y4;
      x4 = x3 + (x2 - x1);
      y4 = f(x4);
      if progflag; print "expand right to " x4 y4; endif;
         goto right;
      endif;
      if progflag; print "contracting"; endif;
         do until dy < tol;
            if (y1 >= y2) and (y1 >= y3) and (y2 <= y3);
               x4 = x2; y4 = y2;
               x2 = x1 + a*(x4 - x1);
               x3 = x4 - a*(x4 - x1);
               y2 = f(x2); y3 = f(x3);
           elseif (y2 > y3) or ((y2 == y3) and (y1 >= y4));
              x4 = x3;  y4 = y3;
              x3 = x2;  y3 = y2;
              x2 = x1 + (x4 - x3);
              y2 = f(x2);
           else;
              x1 = x2;  y1 = y2;
              x2 = x3;  y2 = y3;
              x3 = x4 - (x2 - x1);
              y3 = f(x3);
           endif;
           dy = maxc(y1|y2|y3|y4) - minc(y1|y2|y3|y4);
           if progflag;
              print "x1 x2 x3 x4 dy" x1 x2 x3 x4 dy;
              print "y1 y2 y3 y4   " y1 y2 y3 y4 ;
           endif;
        endo;

        missing = {.};
        if y1 > y2;
           d2y = missing;
           retp (x1, y1, d2y);
        elseif y4 > y3;         @ this case should never occur @
           d2y = missing;
           retp (x4, y4, d2y);
        elseif y2 > y3;
           dx2 = x2-x1;
           dx3 = x3-x1;
           dy2 = y2-y1;
           dy3 = y3-y1;
           d2y = 2*(dx2*dy3-dx3*dy2)/(dx2*dx3*(dx3-dx2));
           retp (x2, y2, d2y);
        else;                  @ y3 >= y2 @
           dx2 = x3-x2;
           dx3 = x4-x2;
           dy2 = y3-y2;
           dy3 = y4-y2;
           d2y = 2*(dx2*dy3-dx3*dy2)/(dx2*dx3*(dx3-dx2));
           retp (x3, y3, d2y);
        endif;
endp;

2 Answers

0

The first step is to get your data from Excel into GAUSS. Here is an example Excel file, example.xls, that you can download for testing.

After you have downloaded the file, you can import the data into GAUSS like this:

x = xlsReadM("example.xls", "B2:B13", 1, "");
y = xlsReadM("example.xls", "C2:C13", 1, "");

If you save the file in a directory other than your GAUSS current working directory, you can add full path to the commands like this:

x = xlsReadM("C:\\mypath\\example.xls", "B2:B13", 1, "");

Inside of quotes, Windows path separators need to be 'escaped', that is why they have two backslashes above. Now you should have the data from the 'B' column of example.xls in x and data from the 'C' column in y.

The next step is to call the function, which should look something like this:

tol = 1e-5;
progflag = 1;

{ bf, pf, bs, ps, u, 
  s2eps, rho, serho, 
  tinf, setinf, logl, lr } = als(x, y, tol, progflag);

After you figure out the data import from Excel, it might be better to try and call the function with your data since some algorithms need more than 12 observations to work properly.

link

aptech

1,773

0

Thank you very much.

link

Tomislav

0