proc quad(xs,h,&u);
@Calculate Quadratic Approximation

     xs=vector containing steady state values.
     h=small positive constant
     u=statement function to be approximated @

local i,j,n,q,z,hz,u:fn;
n=rows(xs);
q=zeros(n+1,n+1);
z=zeros(n,n);
hz=abs(h*xs);
z=diagrv(z,hz);
i=1;
do until i>n;
     if hz[i,1] == 0; z[i,i]=h; endif;
     q[i+1,i+1]=((u(xs+z[.,i])+u(xs-z[.,i]))/2-u(xs))/z[i,i]^2;
     j=1;
     do until j==i;
          q[i+1,j+1]=(u(xs+z[.,i]+z[.,j])-u(xs+z[.,i]-z[.,j])
                     -u(xs-z[.,i]+z[.,j])+u(xs-z[.,i]-z[.,j]))/
                     (8*z[i,i]*z[j,j]);
          q[j+1,i+1]=q[i+1,j+1];
          j=j+1;
     endo;
     i=i+1;
endo;
i=1;
do until i>n;
     q[i+1,1]=(u(xs+z[.,i])-u(xs-z[.,i]))/(2*z[i,i]);
     q[i+1,1]=(q[i+1,1]-2*(q[i+1,2:n+1]*xs))/2;
     q[1,i+1]=q[i+1,1];
     i=i+1;
endo;
retp(q);
endp;

proc reduce(r,a);
@ r is an n x n matrix to be reduced by one row and column
  a is a vector of parameters of a linear equation x(n)=a*[x(1),..,x(n-1)]
      used to eliminate the last row and column of r.  @
local n,x,rnew;
n=rows(r);
x=eye(n-1)|a;
rnew=x'r*x;
retp(rnew);
endp;

proc quadp(xs,h,&u);
@Calculate Quadratic Approximation

     xs=vector containing steady state values.
     h=small positive constant
     u=proc containing function to be approximated @

local i,j,n,q,z,hz,u:proc;
n=rows(xs);
q=zeros(n+1,n+1);
z=zeros(n,n);
hz=abs(h*xs);
z=diagrv(z,hz);
i=1;
do until i>n;
     if hz[i,1] == 0; z[i,i]=h; endif;
     q[i+1,i+1]=((u(xs+z[.,i])+u(xs-z[.,i]))/2-u(xs))/z[i,i]^2;
     j=1;
     do until j==i;
          q[i+1,j+1]=(u(xs+z[.,i]+z[.,j])-u(xs+z[.,i]-z[.,j])
                     -u(xs-z[.,i]+z[.,j])+u(xs-z[.,i]-z[.,j]))/
                     (8*z[i,i]*z[j,j]);
          q[j+1,i+1]=q[i+1,j+1];
          j=j+1;
     endo;
     i=i+1;
endo;
i=1;
do until i>n;
     q[i+1,1]=(u(xs+z[.,i])-u(xs-z[.,i]))/(2*z[i,i]);
     q[i+1,1]=(q[i+1,1]-2*(q[i+1,2:n+1]*xs))/2;
     q[1,i+1]=q[i+1,1];
     i=i+1;
endo;

retp(q);
endp;

proc hpf(n,l);
@    Calculates an n x n matrix that can be used to compute
      cyclical component.
      inputs:   n = number of observations in data set
                l = lambda (see Prescott, "Theory Ahead of Measurement..")
      output:  hpfilt = n x n matrix    @
local nobs,a,i,ai;
nobs=n;
a=zeros(nobs,nobs);
a[1,1]=1+l;
a[1,2]=-2*l; a[1,3]=l;
a[2,1]=-2*l; a[2,2]=1+5*l;
a[2,3]=-4*l; a[2,4]=l;
a[nobs-1,nobs-3]=l;
a[nobs-1,nobs-2]=-4*l;
a[nobs-1,nobs-1]=1+5*l;
a[nobs-1,nobs]=-2*l;
a[nobs,nobs-2]=l;
a[nobs,nobs-1]=-2*l; a[nobs,nobs]=1+l;
i=3;
do until i>(nobs-2);
     a[i,i]=1+6*l;
     a[i,i+1]=-4*l;
     a[i,i+2]=l;
     a[i,i-1]=-4*l;
     a[i,i-2]=l;
     i=i+1;
endo;
ai=eye(nobs)-inv(a);
save \gauss\prog\hpfilt=ai;
retp(ai);
endp;

