0

我正在将演示时态数据分析的 SAS 程序转换为 R。我想使用 R 重现 SAS PROC SPECTRA 输出。

所以，我的问题是，R中的spectrum()函数产生的光谱密度值是否可以转换为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;

到目前为止，在 R 中使用 read.xlsx() 输入数据后，这就是我所得到的：

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

使用 ma 作为 2 阶移动平均平滑器。

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

变量 U 的白噪声测试

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

傅里叶分析的周期图。SAS 产生的周期图值看起来大约是 P（下）乘以 4。这并不奇怪，有人告诉我，某些软件会产生周期图值除以 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")

因此，R 的spectrum() 函数产生的频谱密度值$spec 值与SAS PROC SPECTRA 产生的那些不同。我可以将我的 R 值转换为 SAS 值吗？

就这样。感谢您的时间。

4

0 回答 0