Converting spectral density values produced by spectrum() in R to values produced by SAS PROC SPECTRA

I am converting SAS programs that demonstrate temporal data analysis into R. I would like to reproduce SAS PROC SPECTRA output using R.

So, my question is, can the spectral density values produced by the spectrum() function in R be converted into the values for spectral density produced by SAS PROC SPECTRA?

DATA AR1_09;    
 INPUT t U; 
 OUTPUT;    
CARDS;  
1   -5.19859
2   4.91364
3   -3.86515
4   4.02932
5   -4.12263
6   3.46548
7   -3.01139
8   3.13753
9   -2.34875
10  2.1531
11  -2.01086
12  1.88911
13  -2.22766
14  1.94077
15  0.1786
16  0.84228
17  -1.51301
18  2.62644
19  -3.44148
20  3.13813
21  -2.34959
22  2.70754
23  -2.54789
24  2.04427
25  -2.34041
26  1.13443
27  -0.11853
28  0.74645
29  0.02448
30  0.57811
31  -1.54715
32  1.05646
33  -0.56458
34  0.6863
35  -0.53347
36  0.60813
37  -1.22044
38  0.13136
39  -0.45568
40  0.13459
41  -0.10892
42  0.46324
43  1.01367
44  -2.44015
45  1.62849
46  1.54928
47  -2.7146
48  2.20448
49  -1.58668
50  1.06419
51  -1.41402
52  1.30755
53  -1.55331
54  1.58191
55  -2.38216
56  1.45702
57  0.79562
58  -0.91078
59  -0.59827
60  1.44958
61  -1.81996
62  -0.05101
63  -0.13188
64  1.34861
65  -1.81912
66  0.73641
67  -0.32049
68  -0.37179
69  2.26288
70  -2.2773
71  0.95193
72  -1.24679
73  0.67123
74  -0.40868
75  1.46308
76  -0.71945
77  1.07481
78  -2.25127
79  1.87573
80  -1.52811
81  1.27772
82  -2.96657
83  3.58684
84  -1.7656
85  2.92004
86  -2.36525
87  2.17087
88  -1.65458
89  0.86588
90  0.19505
91  -2.34264
92  3.51124
93  -3.33501
94  3.13522
95  -1.8957
96  0.93527
97  -0.96551
98  0.08307
99  -0.14018
100 0.48641
;   
PROC SPECTRA DATA=AR1_09 OUT=AR1_09PSPEC1 P S WHITETEST;    
 VAR U; 
WEIGHTS 1 2 1;  
RUN;    
PROC SPECTRA DATA=AR1_09 OUT=AR1_09PSPEC2 S WHITETEST;  
 VAR U; 
 WEIGHTS 1 2 3 4 5 4 3 2 1;     
RUN;    
DATA AR1_09PSPEC12; 
 SET AR1_09PSPEC1;  
 n=100; 
 fre=0.5*n*FREQ/(4*ATAN(1));    
 P_01=P_01/(16*ATAN(1));    
 KEEP fre P_01 S_01;    
 RUN;   
DATA AR1_09PSPEC22; 
 SET AR1_09PSPEC2;  
 n=100; 
 fre=0.5*n*FREQ/(4*ATAN(1));    
 S_02=S_01; 
 KEEP fre S_02; 
DATA AR1_09TRUESPEC;    
 SET AR1_09PSPEC1;  
 n=100; 
 rho=-0.9;  
 theoreticalS=1.0/(8*ATAN(1)*(1-2*rho*cos(FREQ)+rho*rho));  
 fre=0.5*n*FREQ/(4*ATAN(1));    
 KEEP fre theoreticalS; 
DATA AR1_09PSPEC;   
 MERGE AR1_09PSPEC12 AR1_09PSPEC22 AR1_09TRUESPEC;  
PROC PRINT DATA=AR1_09PSPEC;    
 VAR fre P_01 S_01 S_02 theoreticalS;   
RUN;    

and so far, after entering the data using read.xlsx() in R, this is what I've got:

AR1.09 <- as.ts(AR1.09[, 2])
install.packages("forecast")
library(forecast)

Using ma for the moving average smoother of order 2.

MA2AR1.09 <- ma(AR1.09, order = 2)
AR1.09PSPEC1 <- spectrum(na.omit(MA2AR1.09))

Tests for White Noise for Variable U

Box.test (MA2AR1.09)
Box.test (MA2AR1.09, type = "Ljung")

Periodogram of Fourier analysis. The periodogram values produced by SAS look to be approximately P (below) times 4. This isn't surprising, I have been told that some software produces periodogram values divided by 4pi.

n <- length(AR1.09)
FF <- abs(fft(AR1.09) / sqrt(n))^2
P <- (4 / n) * FF[1:((n / 2) + 1)]
f <- (0:(n/2)) / n
plot(f, P, type = "h")

So, the $spec values that are the spectral density values produced by R's spectrum() function are not the same as those produced by SAS PROC SPECTRA. Can I transform my R values into the SAS values?

That's all. Thank you for your time.

Upvotes: 0

Views: 118

Answers (0)

Related Questions