Chapter Contents Previous Next
 The SPECTRA Procedure

## Example 17.1: Spectral Analysis of Sunspot Activity

This example analyzes Wolfer's sunspot data (Anderson 1971). The following statements read and plot the data.


title "Wolfer's Sunspot Data";
data sunspot;
input year wolfer @@;
datalines;
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
;

symbol1 i=splines v=dot;
proc gplot data=sunspot;
plot wolfer*year;
run;


The plot of the sunspot series is shown in Output 17.1.1.

Output 17.1.1: Plot of Original Data

The spectral analysis of the sunspot series is performed by the following statements:


proc spectra data=sunspot out=b p s adjmean whitetest;
var wolfer;
weights 1 2 3 4 3 2 1;
run;

proc print data=b(obs=12);
run;


The PROC SPECTRA statement specifies the P and S options to write the periodogram and spectral density estimates to the OUT= data set B. The WEIGHTS statement specifies a triangular spectral window for smoothing the periodogram to produce the spectral density estimate. The ADJMEAN option zeros the frequency 0 value and avoids the need to exclude that observation from the plots. The WHITETEST option prints tests for white noise.

The Fisher's Kappa test statistic of 16.070 is larger than the 5% critical value of 7.2, so the null hypothesis that the sunspot series is white noise is rejected.

The Bartlett's Kolmogorov-Smirnov statistic of 0.6501 is greater than

so reject the null hypothesis that the spectrum represents white noise.

The printed output produced by PROC SPECTRA is shown in Output 17.1.2. The output data set B created by PROC SPECTRA is shown in part in Output 17.1.3.

Output 17.1.2: White Noise Test Results

 Wolfer's Sunspot Data

 SPECTRA Procedure

 Test for White Noise forVariable wolfer M-1 87 Max(P(*)) 4062267 Sum(P(*)) 21156512

 Fisher's Kappa: (M-1)*Max(P(*))/Sum(P(*)) Kappa 16.70489

 Bartlett's Kolmogorov-SmirnovStatistic: Maximum absolute differenceof the standardized partial sums of the periodogramand the CDF of a uniform(0,1) random variable. Test Statistic 0.650055

Output 17.1.3: First 12 Observations of the OUT= Data Set

 Wolfer's Sunspot Data

 Obs FREQ PERIOD P_01 S_01 1 0.00000 . 0.00 59327.52 2 0.03570 176.000 3178.15 61757.98 3 0.07140 88.000 2435433.22 69528.68 4 0.10710 58.667 1077495.76 66087.57 5 0.14280 44.000 491850.36 53352.02 6 0.17850 35.200 2581.12 36678.14 7 0.21420 29.333 181163.15 20604.52 8 0.24990 25.143 283057.60 15132.81 9 0.28560 22.000 188672.97 13265.89 10 0.32130 19.556 122673.94 14953.32 11 0.35700 17.600 58532.93 16402.84 12 0.39270 16.000 213405.16 18562.13

The following statements plot the periodogram and spectral density estimate:


proc gplot data=b;
plot p_01 * freq;
plot p_01 * period;
plot s_01 * freq;
plot s_01 * period;
run;


The periodogram is plotted against frequency in Output 17.1.4 and plotted against period in Output 17.1.5. The spectral density estimate is plotted against frequency in Output 17.1.6 and plotted against period in Output 17.1.7.

Output 17.1.4: Plot of Periodogram by Frequency

Output 17.1.5: Plot of Periodogram by Period

Output 17.1.6: Plot of Spectral Density Estimate by Frequency

Output 17.1.7: Plot of Spectral Density Estimate by Period

Since PERIOD is the reciprocal of frequency, the plot axis for PERIOD is stretched for low frequencies and compressed at high frequencies. One way to correct for this is to use a WHERE statement to restrict the plots and exclude the low frequency components. The following statements plot the spectral density for periods less than 50.


proc gplot data=b;
where period < 50;
plot s_01 * period / HREF=11;
run;


The spectral analysis of the sunspot series confirms a strong 11-year cycle of sunspot activity. The plot makes this clear by drawing a reference line at the 11 year period, which highlights the position of the main peak in the spectral density.

Output 17.1.8 shows the plot. Contrast Output 17.1.8 with Output 17.1.7.

Output 17.1.8: Plot of Spectral Density Estimate by Period to 50 Years

 Chapter Contents Previous Next Top