Article From: Department of Economics
Subject: LGMMQSR.SET

Posted by: Information Provider
Phone Number: (202) 885-3770
E-Mail Address: econ@american.edu

Post Date: 28 Sep 1994
Expiration Date: 28 Oct 2004


               LGMMQSR.SET (Ogaki)

 /* FILE LGMMQSR.SET
      CONDITIONAL HETEROSCEDASTICITY CASE
   Linear restrictions are imposed and tested.
   PROCEDURE LGMM:
   L*m*tend must be less than 8190.
   L>=K shoud hold to be identified.

Model: y(t)=x(t)'b+u(t)
         E[u(t)|I(t)]=0
         z(t) is in I(t).

 Syntax: b=lgmmhetr(y,xp,zp,tend,ru,zero,maxiter,rm,rv,w0flag,st,wav,maxd);

   y=|y(1)   |   y=m*tend by 1 matrix.
     |  |    |   y(t) mx1 matrix, t=1,..,tend.
     |y(tend)|

   xp=|x(1)'   |  xp=m*tend by k matrix.
      |x(2)'   |  x(t) kxm matrix, t=1,..,tend.
      |  |     |
      |x(tend)'|

   zp=|z(1)'   |  zp=m*tend by L matrix.
      |z(2)'   |  z(t) Lxm matrix, t=1,..,tend.
      |  |     |
      |z(tend)'|

   tend     the number of observations
   ru       If u(t) is known to be in I(t+ru+1), give the number ru.
            this is not used if w0flag=2
   w0flag:  a scalar
            if w0flag=1, truncated kernel is used with autocorrelation
                up to the order ru.
            if w0flag=2, QS kernel is used with the following parameters
               to control the kernel.
   st:      a scalar to control the bandwidth parameter for the QS kernel.
         if st=0, then an automatic bandwidth estimator is used.
             if st/=0, then st is used as the bandwidth parameter.
   wav:  (m2+1) by 1 vector; weights given to the a-th element of w(t) for
      the automatic bandwidth estimator (if st=0).
   maxd:    scalar;
             if maxd=0, then nonprewhitened HAC with the QS kernel is used.
         if 0<maxd<1, then the elements of DeltaLS with the absolute
                value greater than maxd is replaced by maxd.  See Andrews
                and Monahan's (1990) footnote 4.
This Program Calculates and Iterates:
 initial condition w=I,
1. sigmazy=1/tend*(sum t=1,tend{z(t)y(t)})
   sigmazx=1/tend*(sum t=1,tend{z(t)x(t)'})
   bu=inv(sigmazx'*w*sigmazx)*sigmazx'w*sigmazy (unrestrictted)
   br:restricted by Rm*b=rv
   u(t)=y(t)-x(t)'br
   if w0flag=1;
       Appropriate elements of Rzu(j) will be replaced by zeros.
    Rzu(j)=1/tend*(sum t=j+1,tend{z(t)u(t)u(t-j)z(t-j)'}  j=0,1,..,ru
    w1=Rzu(0)+{Rzu(1)+Rzu(1)'}+,..,+{Rzu(r)+Rzu(r)'}
    w=inv(w1)
   elseif w0flag=2;
    QS kernel is used, ru is not used in this case.
   endif;
go back to 1.
                */
@ VARQS.PRC, this part defined QS kernel estimates@
proc varqc(wv,np,st,wav,maxd);
local tend,tendm,sigma,wvstr,omegas,tenda,als,va,bls,clsc,dells,indm,i,
      alpha2,alpha2d,iwsws,rhoa,ea,siga4,v2512,v65,omegat,tau,omega,
      xtau,dinv;
   @wv contains w(t) t=1,...,T w(t)=[u(t)' v(t)']' m by 1 vector @
tend=rows(wv);
tendm=tend-np;

 sigma=(wv'wv)/tendm;
if maxd==0;
 ?"Nonprewhitened HAC with the QS kernel is used for calculation of omega";
 wvstr=wv;
 omegas=sigma;
 tenda=tend;
else;
 ? "Prewhitened HAC with the QS kernel is used for calculation of omega";
 tenda=tend-1;
 als=wv[2:tend,.]'wv[1:tend-1,.]*invpd(wv[1:tend-1,.]'wv[1:tend-1,.]);
 {va,bls}=eigrs2(als*(als'));
 {va,clsc}=eigrs2(als'als);
 dells=bls'als*clsc;
"ALS and DELTALS for prewhitening=" als dells;
"maxd=" maxd;
 indm=0;
 i=1;
 do until i>rows(als);
  if dells[i,i]>maxd;
   ? i "-th element of DeltaLS=" dells[i,i] "is replaced by" maxd;
   indm=1;
   dells[i,i]=maxd;
  elseif dells[i,i]<-maxd;
   ? i "-th element of DeltaLS=" dells[i,i] "is replaced by" -maxd;
   indm=1;
   dells[i,i]=-maxd;
  endif;
  i=i+1;
 endo;
 if indm==1;
  als=bls*dells*(clsc');
 endif;
 wvstr=wv[2:tend,.]-wv[1:tend-1,.]*(als');
       @ a bug found on 4/30/91, used to be "*als" @
 omegas=(wvstr'wvstr)/tendm;
endif;
if st==0;
? "***** Automatic Bandwidth Estimator is used";
 clear alpha2,alpha2d;
 i=1;
 do until i>cols(wv);
  if wav[i,1]>0;
   iwsws=invpd(wvstr[1:tenda-1,i]'wvstr[1:tenda-1,i]);
   rhoa=iwsws*(wvstr[1:tenda-1,i]'wvstr[2:tenda,i]);
   ea=wvstr[2:tenda,i]-wvstr[1:tenda-1,i]*rhoa;
   siga4=((ea'ea)/(tenda-1))^2;
   alpha2=alpha2+4*wav[i,1]*(rhoa^2)*(siga4)/((1-rhoa)^8);
   alpha2d=alpha2d+wav[i,1]*siga4/((1-rhoa)^4);
  endif;
  i=i+1;
 endo;
 alpha2=alpha2/alpha2d;
 st=1.3221*((alpha2*(tend))^0.2);
endif;
 ? "The bandwidth parameter used" st;
 v2512=25/(12*(pi^2));
 v65=6*pi/5;
 tau=1;
 do until tau>tenda-1;           @change from >= 5/14/91@
  omegat=(wvstr[1:tenda-tau,.]'wvstr[1+tau:tenda,.])/tendm;
  xtau=tau/st;
  omegas=omegas+v2512*(sin(v65*xtau)/(v65*xtau)-cos(v65*xtau))
                     *(omegat+omegat')/(xtau^2);
  tau=tau+1;
 endo;
if maxd==0;
 omega=omegas;
else;
 dinv=inv(eye(rows(als))-als);
 omega=dinv*omegas*(dinv');
endif;
retp(omega);
endp;

proc lgmmhetr(y,xp,zp,tend,ru,zero,maxiter,rm,rv,w0flag,st,wav,maxd);
local sigmazy,sigmazx,differ,iter,u,v,se,g,prob,t,ruj,rzuj,vt,ut,
      zpt,zptj,w,j,i,vtj,utj,df,chi,w0,br,bu,varbr,
      varbu,zxw,tstr,rpri,rbur,k,l,m,wv;
clear chi;
@"======== Linear Restrictions Rm*b=rv";
 ? "Rm=" rm;
 ? "rv'=" rv';@
k=cols(xp);
L=cols(zp);
m=round(rows(y)/tend);
sigmazy=zp'y/tend;
sigmazx=zp'xp/tend;
w0=eye(l);
df=L-k;
differ=10;iter=1;
do until (differ<zero) or (iter>maxiter);
 zxw=sigmazx'w0;
 v=invpd(zxw*sigmazx);
 bu=v*zxw*sigmazy;   @unrestricted@

 rpri=invpd(rm*v*(rm'));
 rbur=rm*bu-rv;
 br=bu-v*(rm')*rpri*rbur;

 u=y-xp*br;

 varbu=v/tend;
 varbr=(v-v*(rm')*rpri*rm*v)/tend;
 tstr=(rbur')*invpd(rm*varbu*(rm'))*rbur;
 save br,varbr,bu;

if iter==1;
 goto lab10;
endif;
 se=sqrt(diag(varbr));
 g=zp'u/tend;
 if df>0;
 chi=tend*(g'w0*g);
 prob=cdfchic(chi,df);endif;
 ? "********************* ITERATION=" iter "********************";
 ? "restricted  br=" br';
 ? "          s.e.=" se';
 ? "unrestricted  bu=" bu';
 ? "            s.e.=" sqrt(diag(varbu))';
 ? "LM Test for Rm*b=rv";
 tstr "(" cdfchic(tstr,rows(rm)) ") df=" rows(rm);

 if df>0;
 ? "   chi square for overidentifying restrictions=";
 chi "(" prob ")";
 elseif df==0;? "   Just Identified";
 else;? "   L<k; Not Identified";endif;
 ? "   d.f.=" df;
lab10:

if w0flag==1;
? "******  Truncated Kernel Used  *****";
 t=1;clear rzuj;
 do until t>tend;
  vt=seqa(m*(t-1)+1,1,m);
  zpt=zp[vt,0];ut=u[vt,0];
  rzuj=rzuj+(zpt')*ut*(ut')*zpt;
  t=t+1;
 endo;
 w=rzuj;
 j=1;
 do until j>ru;
  t=j+1;clear rzuj;
  do until t>tend;
   vt=seqa(m*(t-1)+1,1,m);vtj=seqa(m*(t-j-1)+1,1,m);
   ut=u[vt,0];utj=u[vtj,0];
   zpt=zp[vt,0];zptj=zp[vtj,0];
   rzuj=rzuj+(zpt')*ut*(utj')*zptj;
   t=t+1;
  endo;
  w=w+rzuj+(rzuj');
  j=j+1;
 endo;
 w=invpd(w/tend);
 differ=w-w0;differ=maxc(maxc(abs(differ)));
 w0=w;

elseif w0flag==2;
   @ **** QS kernel ****@
 t=1;
 wv=zeros(tend,L);
 do until t>tend;
  vt=seqa(m*(t-1)+1,1,m);
  zpt=zp[vt,0];ut=u[vt,0];
  wv[t,.]=(ut')*zpt;
  t=t+1;
 endo;
 w=invpd(varqc(wv,k,st,wav,maxd));
 differ=w-w0;differ=maxc(maxc(abs(differ)));
 w0=w;
endif;

 iter=iter+1;
 ? "max(|differ|)=" differ;
endo;

retp(chi);
endp;

************************************************************************
Department of Economics is entirely responsible for the information
provided above. Please direct any comments to Information Provider at:
          E-Mail Address: econ@american.edu
          Phone Number:   (202) 885-3770