/* Appendix 1: Calculating P-Values for Fisher's Kappa */
/* Calculate p-values for various sample */
/* sizes and values of Fisher's Kappa. */
data one;
TOOBIG=1e50; /* Machine upper limit on floating point */
TOOSMALL=1e-50; /* Machine lower limit on floating point */
MAXFACT=41; /* Maximum factorial for GAMMA function */
retain sum;
/* n is the series length */
do n=25,50,100,200,500,1000,3000;
m=floor(n/2); /* 1-based indexing uses n */
/* rather than n-1 */
/* r is the rank of the periodogram ordinate */
do r=1 to 4;
/* z represents a possible value of Fisher's */
/* Kappa. Modify this loop to include values */
/* of interest. */
do z=0.05 to 0.95 by 0.05;
b=floor(1/z);
if b<(1/z) then a=b;
else a=b-1;
gprob=0;
top=min(a,m); /* Largest index for rank in */
/* summation */
sum=1.0;
do j=r to top;
if j=r then do;
/* Initialize for first step */
dflag=0;
/* Use GAMMA function rather than */
/* multiplication if possible. */
if r<MAXFACT then do;
denom=gamma(r)*r;
prob=m/denom;
end;
else do;
/* dflag=1 implies that we must now keep */
/* track of denominator. */
dflag=1;
denom1=max(r-1,1);
denom=r*denom1;
prob=m;
end;
nmult=m-1;
amult=1-r*z;
bot1=m-r; /* Lower limit of mult. */
curexp=m-1; /* Current exponent value */
do while(nmult>bot1);
bot2=TOOBIG/nmult;
/* Loop until overflow avoided */
do while(prob>bot2);
/* If it is impossible to avoid */
/* overflow, set prob to missing. */
/* Use -1 for now as a flag, */
/* since comparisons using */
/* missing values can be risky. */
if curexp=0 and dflag=0
then prob=-1;
/* Instead of raising to a power, */
/* multiply one term at a time */
/* and anticipate using a */
/* multiplication or division to */
/* prevent overflow or underflow. */
else if curexp>0 then do;
prob=prob*amult; /* Cannot */
/* overflow */
curexp=curexp-1;
end;
else if dflag=1 and
prob>denom*TOOSMALL
then do;
prob=prob/denom; /* Cannot */
/* underflow */
if denom1>1 then
denom1=denom1-1;
if denom1=1 then
dflag=0;
denom=denom1;
end;
end;
/* Failure */
if prob=-1 then do;
prob=.;
nmult=0;
dflag=0;
end;
else do;
/* Cannot overflow */
prob=prob*nmult;
nmult=nmult-1;
end;
end;
/* Now it is okay to exponentiate */
if curexp>0 then
prob=prob*(amult**curexp);
if dflag=1 then
do while(dflag=1);
prob=prob/denom;
if denom1>1 then
denom1=denom1-1;
if denom1=1 then
dflag=0;
denom=denom1;
end;
sum=prob;
end;
else do; /* Subsequent steps in sum for j>r */
/* using recursion. We are less */
/* careful because overflow is */
/* unlikely and underflow is well */
/* approximated by zero. */
f1=(1-j*z)/(1-(j-1)*z);
f2=(j-1)*(m-j+1)/(j*(j-r));
prob=f2*(f1**(m-1));
sum=-sum*prob;
end;
gprob=gprob+sum;
end; /* End j loop */
output;
end; /* End z loop */
end; /* End r loop */
end; /* End n loop */
keep z n m r gprob;
run;
proc print data=one;
run;
/**********************************************************************/
/* Appendix 2: Calculating P-Values for Fisher's Kappa */
/* from PROC SPECTRA Output */
/* Data from SAS/ETS Sample Library: SPECTRA1 */
data sunspot;
input year wolfer @@;
cards;
1749 809 1750 834 1751 477 1752 478 1753 307 1754 122 1755 96
1756 102 1757 324 1758 476 1759 540 1760 629 1761 859 1762 612
1763 451 1764 364 1765 209 1766 114 1767 378 1768 698 1769 1061
1770 1008 1771 816 1772 665 1773 348 1774 306 1775 70 1776 198
1777 925 1778 1544 1779 1259 1780 848 1781 681 1782 385 1783 228
1784 102 1785 241 1786 829 1787 1320 1788 1309 1789 1181 1790 899
1791 666 1792 600 1793 469 1794 410 1795 213 1796 160 1797 64
1798 41 1799 68 1800 145 1801 340 1802 450 1803 431 1804 475
1805 422 1806 281 1807 101 1808 81 1809 25 1810 0 1811 14
1812 50 1813 122 1814 139 1815 354 1816 458 1817 411 1818 304
1819 239 1820 157 1821 66 1822 40 1823 18 1824 85 1825 166
1826 363 1827 497 1828 625 1829 670 1830 710 1831 478 1832 275
1833 85 1834 132 1835 569 1836 1215 1837 1383 1838 1032 1839 858
1840 632 1841 368 1842 242 1843 107 1844 150 1845 401 1846 615
1847 985 1848 1243 1849 959 1850 665 1851 645 1852 542 1853 390
1854 206 1855 67 1856 43 1857 228 1858 548 1859 938 1860 957
1861 772 1862 591 1863 440 1864 470 1865 305 1866 163 1867 73
1868 373 1869 739 1870 1391 1871 1112 1872 1017 1873 663 1874 447
1875 171 1876 113 1877 123 1878 34 1879 60 1880 323 1881 543
1882 597 1883 637 1884 635 1885 522 1886 254 1887 131 1888 68
1889 63 1890 71 1891 356 1892 730 1893 849 1894 780 1895 640
1896 418 1897 262 1898 267 1899 121 1900 95 1901 27 1902 50
1903 244 1904 420 1905 635 1906 538 1907 620 1908 485 1909 439
1910 186 1911 57 1912 36 1913 14 1914 96 1915 474 1916 571
1917 1039 1918 806 1919 636 1920 376 1921 261 1922 142 1923 58
1924 167
;
proc spectra data=sunspot p out=a whitetest;
var wolfer;
run;
data a;
set a;
retain n;
if period=. then delete;
if period=2.0 then delete;
if _n_=2 then n=period;
run;
proc means data=a sum noprint;
var p_01;
id n;
output out=b sum=psum;
run;
proc sort data=a;
by descending p_01;
run;
data c;
merge a b;
by n;
drop _type_ _freq_;
run;
data _null_;
set c;
TOOBIG=1e50; /* Machine upper limit on floating point */
TOOSMALL=1e-50; /* Machine lower limit on floating point */
MAXFACT=41; /* Maximum Factorial */
m=floor((n-1)/2);
r=_n_;
if r>12 then delete; /* Upper limit on r */
z=p_01/psum;
b=floor(1/z);
if b<(1/z) then a=b;
else a=b-1;
gprob=0;
top=min(a,m);
do j=r to top;
if j=r then do; /* Initialize for first step */
dflag=0;
/* Use GAMMA function rather than multiplication */
/* if possible. */
if r<MAXFACT then do;
denom=gamma(r)*r;
prob=m/denom;
end;
else do;
/* dflag=1 implies that we must now keep */
/* track of denominator. */
dflag=1;
denom1=max(r-1,1);
denom=r*denom1;
prob=m;
end;
nmult=m-1;
amult=1-r*z;
bot1=m-r; /* Lower limit of mult. */
curexp=m-1; /* Current exponent value */
do while(nmult>bot1);
bot2=TOOBIG/nmult;
/* Loop until overflow avoided */
do while(prob>bot2);
/* If it is impossible to avoid overflow, */
/* set prob to missing. Use -1 for now as */
/* a flag, since comparisons using missing */
/* values can be risky. */
if curexp=0 and dflag=0
then prob=-1;
/* Instead of raising to a power, multiply */
/* one term at a time and anticipate using */
/* a multiplication or division to prevent */
/* overflow or underflow. */
else if curexp>0 then do;
/* Cannot overflow */
prob=prob*amult;
curexp=curexp-1;
end;
else if dflag=1 and
prob>denom*TOOSMALL then do;
/* Cannot underflow */
prob=prob/denom;
if denom1>1 then
denom1=denom1-1;
if denom1=1 then dflag=0;
denom=denom1;
end;
end;
if prob=-1 then do; /* Failure */
prob=.;
nmult=0;
dflag=0;
end;
else do;
/* Cannot overflow */
prob=prob*nmult;
nmult=nmult-1;
end;
end;
/* Now it is okay to exponentiate */
if curexp>0 then
prob=prob*(amult**curexp);
if dflag=1 then do while(dflag=1);
prob=prob/denom;
if denom1>1 then denom1=denom1-1;
if denom1=1 then dflag=0;
denom=denom1;
end;
sum=prob;
end;
else do; /* Subsequent steps in sum for j>r using */
/* recursion. We are less careful because */
/* overflow is unlikely, and underflow is */
/* well approximated by zero. */
f1=(1-j*z)/(1-(j-1)*z);
f2=(j-1)*(m-j+1)/(j*(j-r));
prob=f2*(f1**(m-1));
sum=-sum*prob;
end;
gprob=gprob+sum;
end; /* End j loop */
/* Flag floating-point error */
if gprob>1 then gprob=2.0;
/* Flag floating-point error */
if gprob<0 then gprob=-1.0;
mz=z*m;
if _n_=1 then
put @1 'PERIOD R Z KAPPA P-Value';
put @1 period 6.3 @8 r 3.0 @12 z 10.5
@24 mz 10.5 @36 gprob 8.5;
run;