Computing Roots of Polynomials
/****************************************************************/
/* S A S S A M P L E L I B R A R Y */
/* */
/* NAME: POLYROOT */
/* TITLE: Computing Roots of Polynomials */
/* PRODUCT: IML */
/* SYSTEM: ALL */
/* KEYS: MATH */
/* PROCS: IML */
/* DATA: Some polynomials selected from */
/* */
/* Woodfield, Terry J. (1988). Simulating Stationary Gaussian */
/* ARMA Time Series. Computer Science and Statistics: */
/* Proceedings of the 20th Symposium on the Interface, */
/* Reston, Virginia. */
/* */
/* SUPPORT: TJW UPDATE: 18Sep89 */
/* REF: */
/* */
/****************************************************************/
options ls=80 ps=60;
proc iml;
reset fuzz noname fw=8;
/*--------------------------------------------*/
/*---- Multiply two complex polynomials ----*/
/*--------------------------------------------*/
start prodpoly(ain,bin,c);
k=nrow(ain);
m=nrow(bin);
if ncol(ain)^=2 then do;
if ncol(ain)=1 then a=ain||j(k,1,0);
else do;
print 'ERROR: Factor matrix must be N by 1 or N by 2.';
return;
end;
end;
else a=ain;
if ncol(bin)^=2 then do;
if ncol(bin)=1 then b=bin||j(m,1,0);
else do;
print 'ERROR: Factor matrix must be N by 1 or N by 2.';
return;
end;
end;
else b=bin;
n=k+m-1;
c=j(n,2,0);
do i=1 to n;
sumr=0;
sumi=0;
mn=max(1,i-m+1);
mx=min(i,k);
do j=mn to mx;
sumr=sumr+a[j,1]*b[i+1-j,1]-a[j,2]*b[i+1-j,2];
sumi=sumi+a[j,2]*b[i+1-j,1]+a[j,1]*b[i+1-j,2];
end;
c[i,1]=sumr;
c[i,2]=sumi;
end;
free i j k m n mn mx sumr sumi a b;
finish; /*---- End of prodpoly ----*/
/*--------------------------------------*/
/*---- Form polynomial from roots ----*/
/*--------------------------------------*/
start formpoly(roots,poly,tol);
if nrow(tol)=0 then tol=1e-10;
if ncol(roots)^=2 then do;
if ncol(roots)=1 then rts=roots||j(nrow(roots),1,0);
else do;
print 'ERROR: ROOTS matrix must be N by 2.';
return;
end;
end;
else rts=roots;
n=nrow(rts);
poly1=j(2,2,0);
poly2=j(2,2,0);
poly2[1,1]=1;
poly2[1,2]=0;
poly2[2,1]=-rts[1,1];
poly2[2,2]=-rts[1,2];
do i=2 to n;
poly1[1,1]=1;
poly1[1,2]=0;
poly1[2,1]=-rts[i,1];
poly1[2,2]=-rts[i,2];
run prodpoly(poly1,poly2,c);
poly2=c;
end;
n=n+1;
do i=1 to n;
if abs(poly2[i,2])>tol then do;
print
'ERROR: Complex ROOTS must come in conjugate pairs.';
return;
end;
end;
poly=(poly2[,1])`;
free i n poly1 poly2 tol rts;
finish; /*---- End of formpoly ----*/
/*-------------------------------*/
/*---- Evaluate polynomial ----*/
/*-------------------------------*/
start evalpl(poly,arg,result);
n=max(nrow(poly),ncol(poly));
nargs=nrow(arg);
result=j(nargs,2,0);
do i=1 to nargs;
real=arg[i,1];
imag=arg[i,2];
preal=poly[n];
pimag=0.;
preal=preal+real*poly[n-1] ;
pimag=pimag+imag*poly[n-1] ;
do j=n-2 to 1 by -1;
creal=real;
real=real*arg[i,1]-imag*arg[i,2] ;
imag=imag*arg[i,1]+creal*arg[i,2] ;
preal=preal+real*poly[j] ;
pimag=pimag+imag*poly[j] ;
end;
result[i,1]=preal;
result[i,2]=pimag;
end;
print 'Results of evaluating polynomial with coefficients';
if nrow(poly)>ncol(poly) then print (poly`);
else print poly;
tabl=arg||result||j(nargs,2,0);
do i=1 to nargs;
tabl[i,5]=arg[i,1]*arg[i,1]+arg[i,2]*arg[i,2];
tabl[i,6]=result[i,1]*result[i,1]+result[i,2]*result[i,2];
end;
cname={'R(ARG)' 'I(ARG)' 'R(P(ARG))' 'I(P(ARG))'
'|ARG|**2' '|P(ARG)|**2'};
print tabl [colname=cname];
free i j real imag preal pimag creal tabl cname;
finish; /*---- End of evalpl ----*/
/*------------------------------------------*/
/*---- Reverse order of vector values ----*/
/*------------------------------------------*/
start flip(a,b);
n=ncol(a);
b=a[-(-n:-1)];
if nrow(b)>1 then b=b`;
free n;
finish; /*---- End of flip ----*/
/*---------------------------------*/
/*---- MAIN Calculation area ----*/
/*---------------------------------*/
/*---- Roots: 0.9, -1.3 ----*/
roots={0.9 0, -1.3 0};
run formpoly(roots,phi,tol);
b=polyroot(phi);
run evalpl(phi,b,result);
/*---- Roots: 2.4 +/- 5i ----*/
roots={2.4 5, 2.4 -5};
run formpoly(roots,theta,tol);
theta=3*theta;
b=polyroot(theta);
run evalpl(theta,b,result);
/*---- Roots: 0.9, -1.3, 2.4 +/- 5i ----*/
roots={0.9 0, -1.3 0, 2.4 5, 2.4 -5};
run formpoly(roots,phi,tol);
phi=3*phi;
b=polyroot(phi);
run evalpl(phi,b,result);
/*---------------------------------------------------------*/
/*---- The following six models are taken from the ----*/
/*---- Woodfield(1988) paper, "Simulating Univariate ----*/
/*---- Gaussian ARMA Time Series". ----*/
/*---------------------------------------------------------*/
/*---- Model 1 ----*/
phi={1 -0.8 0.65};
run flip(phi,a);
b=polyroot(a);
run evalpl(a,b,r);
/*---- Model 2 ----*/
theta={1 1.25 0.35};
run flip(theta,a);
b=polyroot(a);
run evalpl(a,b,r);
/*---- Model 3 ----*/
phi={1 -0.95};
run flip(phi,a);
b=polyroot(a);
run evalpl(a,b,r);
theta={1 0.85};
run flip(theta,a);
b=polyroot(a);
run evalpl(a,b,r);
/*---- Model 4 ----*/
phi={1 -1.5 1.21 -0.455};
run flip(phi,a);
b=polyroot(a);
run evalpl(a,b,r);
theta={1 0.2 0.9};
run flip(theta,a);
b=polyroot(a);
run evalpl(a,b,r);
/*---- Model 5 ----*/
phi={1 -1 0.24};
run flip(phi,a);
b=polyroot(a);
run evalpl(a,b,r);
theta={1 0.4 0.2 0.1};
run flip(theta,a);
b=polyroot(a);
run evalpl(a,b,r);
/*---- Model 6 ----*/
phi={1 -0.3357 0.0821 0.1570 0.2567};
run flip(phi,a);
b=polyroot(a);
run evalpl(a,b,r);
theta={1 -0.6077 0.0831 0.1903};
run flip(theta,a);
b=polyroot(a);
run evalpl(a,b,r);
quit;
/*-------------------------------*/
/*---- End of POLYROOT.SAS ----*/
/*-------------------------------*/