R version 2.15.2 (2012-10-26) -- "Trick or Treat"
Copyright (C) 2012 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: i686-pc-linux-gnu (32-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> x <- array(list(7645
+ ,27276
+ ,111528
+ ,27383
+ ,7240
+ ,27310
+ ,111627
+ ,27474
+ ,7237
+ ,27346
+ ,111733
+ ,27610
+ ,7170
+ ,27385
+ ,111849
+ ,27460
+ ,7067
+ ,27422
+ ,111970
+ ,27079
+ ,7149
+ ,27459
+ ,112093
+ ,26726
+ ,6979
+ ,27498
+ ,112222
+ ,25414
+ ,6766
+ ,27541
+ ,112354
+ ,25914
+ ,6850
+ ,27584
+ ,112486
+ ,26305
+ ,6731
+ ,27627
+ ,112493
+ ,25660
+ ,6927
+ ,27640
+ ,112596
+ ,26115
+ ,7116
+ ,27666
+ ,112619
+ ,26458
+ ,11299
+ ,27675
+ ,112695
+ ,26722
+ ,10544
+ ,27704
+ ,112737
+ ,26134
+ ,10083
+ ,27709
+ ,112803
+ ,26203
+ ,9501
+ ,27746
+ ,112852
+ ,26568
+ ,9450
+ ,27780
+ ,112912
+ ,26718
+ ,8950
+ ,27816
+ ,113029
+ ,26560
+ ,8578
+ ,27854
+ ,113154
+ ,25473
+ ,8395
+ ,27896
+ ,113281
+ ,25676
+ ,7631
+ ,27939
+ ,113414
+ ,25859
+ ,7816
+ ,27982
+ ,113546
+ ,25632
+ ,7491
+ ,28021
+ ,113573
+ ,25824
+ ,7678
+ ,28052
+ ,113660
+ ,26097
+ ,15124
+ ,28059
+ ,113666
+ ,26179
+ ,15227
+ ,28085
+ ,113758
+ ,25788
+ ,15421
+ ,28118
+ ,113769
+ ,26277
+ ,15012
+ ,28153
+ ,113857
+ ,26342
+ ,14861
+ ,28184
+ ,113953
+ ,26592
+ ,14646
+ ,28217
+ ,114060
+ ,26679
+ ,14727
+ ,28252
+ ,114173
+ ,25694
+ ,14505
+ ,28290
+ ,114288
+ ,26024
+ ,13796
+ ,28330
+ ,114411
+ ,26038
+ ,13389
+ ,28369
+ ,114530
+ ,25736
+ ,12860
+ ,28404
+ ,114632
+ ,25930
+ ,12049
+ ,28437
+ ,114648
+ ,26265
+ ,14393
+ ,28526
+ ,114728
+ ,26042
+ ,15104
+ ,28559
+ ,114735
+ ,24925
+ ,14636
+ ,28591
+ ,114821
+ ,25552
+ ,14574
+ ,28624
+ ,114910
+ ,26142
+ ,14735
+ ,28653
+ ,115001
+ ,26474
+ ,14609
+ ,28685
+ ,115102
+ ,26630
+ ,14517
+ ,28718
+ ,115207
+ ,25460
+ ,14876
+ ,28755
+ ,115317
+ ,25472
+ ,15221
+ ,28794
+ ,115433
+ ,25325
+ ,15128
+ ,28831
+ ,115542
+ ,25121
+ ,15039
+ ,28865
+ ,115640
+ ,25313
+ ,14953
+ ,28896
+ ,115731
+ ,25535
+ ,13097
+ ,28947
+ ,115828
+ ,25255
+ ,13323
+ ,28976
+ ,115907
+ ,24856
+ ,13759
+ ,29005
+ ,115988
+ ,25289
+ ,13897
+ ,29035
+ ,116067
+ ,25411
+ ,13920
+ ,29063
+ ,116156
+ ,25372
+ ,13908
+ ,29093
+ ,116250
+ ,25322
+ ,14024
+ ,29123
+ ,116347
+ ,24961
+ ,13892
+ ,29158
+ ,116453
+ ,24973
+ ,13792
+ ,29193
+ ,116559
+ ,25241
+ ,13628
+ ,29228
+ ,116664
+ ,24832
+ ,13751
+ ,29259
+ ,116755
+ ,24927
+ ,13919
+ ,29286
+ ,116808
+ ,25024
+ ,12258
+ ,29727
+ ,116832
+ ,25131
+ ,12088
+ ,29760
+ ,116896
+ ,24657
+ ,12544
+ ,29792
+ ,116986
+ ,24835
+ ,12794
+ ,29824
+ ,117081
+ ,25126
+ ,12749
+ ,29854
+ ,117177
+ ,25481
+ ,12720
+ ,29885
+ ,117277
+ ,25321
+ ,12500
+ ,29918
+ ,117381
+ ,24787
+ ,12673
+ ,29954
+ ,117492
+ ,24626
+ ,12806
+ ,29991
+ ,117600
+ ,24850
+ ,12758
+ ,30027
+ ,117710
+ ,24566)
+ ,dim=c(4
+ ,70)
+ ,dimnames=list(c('unemployment'
+ ,'black'
+ ,'males'
+ ,'highschool
')
+ ,1:70))
> y <- array(NA,dim=c(4,70),dimnames=list(c('unemployment','black','males','highschool
'),1:70))
> for (i in 1:dim(x)[1])
+ {
+ for (j in 1:dim(x)[2])
+ {
+ y[i,j] <- as.numeric(x[i,j])
+ }
+ }
> par3 = 'No Linear Trend'
> par2 = 'Do not include Seasonal Dummies'
> par1 = '1'
> library(lattice)
> library(lmtest)
Loading required package: zoo
Attaching package: 'zoo'
The following object(s) are masked from 'package:base':
as.Date, as.Date.numeric
> n25 <- 25 #minimum number of obs. for Goldfeld-Quandt test
> par1 <- as.numeric(par1)
> x <- t(y)
> k <- length(x[1,])
> n <- length(x[,1])
> x1 <- cbind(x[,par1], x[,1:k!=par1])
> mycolnames <- c(colnames(x)[par1], colnames(x)[1:k!=par1])
> colnames(x1) <- mycolnames #colnames(x)[par1]
> x <- x1
> if (par3 == 'First Differences'){
+ x2 <- array(0, dim=c(n-1,k), dimnames=list(1:(n-1), paste('(1-B)',colnames(x),sep='')))
+ for (i in 1:n-1) {
+ for (j in 1:k) {
+ x2[i,j] <- x[i+1,j] - x[i,j]
+ }
+ }
+ x <- x2
+ }
> if (par2 == 'Include Monthly Dummies'){
+ x2 <- array(0, dim=c(n,11), dimnames=list(1:n, paste('M', seq(1:11), sep ='')))
+ for (i in 1:11){
+ x2[seq(i,n,12),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> if (par2 == 'Include Quarterly Dummies'){
+ x2 <- array(0, dim=c(n,3), dimnames=list(1:n, paste('Q', seq(1:3), sep ='')))
+ for (i in 1:3){
+ x2[seq(i,n,4),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> k <- length(x[1,])
> if (par3 == 'Linear Trend'){
+ x <- cbind(x, c(1:n))
+ colnames(x)[k+1] <- 't'
+ }
> x
unemployment black males highschool\r
1 7645 27276 111528 27383
2 7240 27310 111627 27474
3 7237 27346 111733 27610
4 7170 27385 111849 27460
5 7067 27422 111970 27079
6 7149 27459 112093 26726
7 6979 27498 112222 25414
8 6766 27541 112354 25914
9 6850 27584 112486 26305
10 6731 27627 112493 25660
11 6927 27640 112596 26115
12 7116 27666 112619 26458
13 11299 27675 112695 26722
14 10544 27704 112737 26134
15 10083 27709 112803 26203
16 9501 27746 112852 26568
17 9450 27780 112912 26718
18 8950 27816 113029 26560
19 8578 27854 113154 25473
20 8395 27896 113281 25676
21 7631 27939 113414 25859
22 7816 27982 113546 25632
23 7491 28021 113573 25824
24 7678 28052 113660 26097
25 15124 28059 113666 26179
26 15227 28085 113758 25788
27 15421 28118 113769 26277
28 15012 28153 113857 26342
29 14861 28184 113953 26592
30 14646 28217 114060 26679
31 14727 28252 114173 25694
32 14505 28290 114288 26024
33 13796 28330 114411 26038
34 13389 28369 114530 25736
35 12860 28404 114632 25930
36 12049 28437 114648 26265
37 14393 28526 114728 26042
38 15104 28559 114735 24925
39 14636 28591 114821 25552
40 14574 28624 114910 26142
41 14735 28653 115001 26474
42 14609 28685 115102 26630
43 14517 28718 115207 25460
44 14876 28755 115317 25472
45 15221 28794 115433 25325
46 15128 28831 115542 25121
47 15039 28865 115640 25313
48 14953 28896 115731 25535
49 13097 28947 115828 25255
50 13323 28976 115907 24856
51 13759 29005 115988 25289
52 13897 29035 116067 25411
53 13920 29063 116156 25372
54 13908 29093 116250 25322
55 14024 29123 116347 24961
56 13892 29158 116453 24973
57 13792 29193 116559 25241
58 13628 29228 116664 24832
59 13751 29259 116755 24927
60 13919 29286 116808 25024
61 12258 29727 116832 25131
62 12088 29760 116896 24657
63 12544 29792 116986 24835
64 12794 29824 117081 25126
65 12749 29854 117177 25481
66 12720 29885 117277 25321
67 12500 29918 117381 24787
68 12673 29954 117492 24626
69 12806 29991 117600 24850
70 12758 30027 117710 24566
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) black males `highschool\\r`
-4.490e+05 -1.009e+01 6.162e+00 1.642e+00
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-3387.3 -1321.7 -418.8 1327.9 4456.6
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -4.490e+05 6.260e+04 -7.172 8.01e-10 ***
black -1.009e+01 2.046e+00 -4.933 5.77e-06 ***
males 6.162e+00 9.529e-01 6.466 1.42e-08 ***
`highschool\\r` 1.642e+00 5.448e-01 3.013 0.00367 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1896 on 66 degrees of freedom
Multiple R-squared: 0.6252, Adjusted R-squared: 0.6082
F-statistic: 36.7 on 3 and 66 DF, p-value: 4.493e-14
> if (n > n25) {
+ kp3 <- k + 3
+ nmkm3 <- n - k - 3
+ gqarr <- array(NA, dim=c(nmkm3-kp3+1,3))
+ numgqtests <- 0
+ numsignificant1 <- 0
+ numsignificant5 <- 0
+ numsignificant10 <- 0
+ for (mypoint in kp3:nmkm3) {
+ j <- 0
+ numgqtests <- numgqtests + 1
+ for (myalt in c('greater', 'two.sided', 'less')) {
+ j <- j + 1
+ gqarr[mypoint-kp3+1,j] <- gqtest(mylm, point=mypoint, alternative=myalt)$p.value
+ }
+ if (gqarr[mypoint-kp3+1,2] < 0.01) numsignificant1 <- numsignificant1 + 1
+ if (gqarr[mypoint-kp3+1,2] < 0.05) numsignificant5 <- numsignificant5 + 1
+ if (gqarr[mypoint-kp3+1,2] < 0.10) numsignificant10 <- numsignificant10 + 1
+ }
+ gqarr
+ }
[,1] [,2] [,3]
[1,] 3.221244e-04 6.442487e-04 9.996779e-01
[2,] 3.607186e-05 7.214372e-05 9.999639e-01
[3,] 2.445644e-06 4.891287e-06 9.999976e-01
[4,] 3.470107e-07 6.940213e-07 9.999997e-01
[5,] 9.689314e-08 1.937863e-07 9.999999e-01
[6,] 4.079102e-08 8.158205e-08 1.000000e+00
[7,] 4.978421e-02 9.956841e-02 9.502158e-01
[8,] 8.461928e-02 1.692386e-01 9.153807e-01
[9,] 6.406989e-02 1.281398e-01 9.359301e-01
[10,] 4.106998e-02 8.213997e-02 9.589300e-01
[11,] 2.913015e-02 5.826031e-02 9.708698e-01
[12,] 2.695256e-02 5.390512e-02 9.730474e-01
[13,] 1.687856e-02 3.375712e-02 9.831214e-01
[14,] 1.552144e-02 3.104289e-02 9.844786e-01
[15,] 4.756650e-02 9.513301e-02 9.524335e-01
[16,] 1.145525e-01 2.291050e-01 8.854475e-01
[17,] 4.752999e-01 9.505998e-01 5.247001e-01
[18,] 9.993189e-01 1.362282e-03 6.811412e-04
[19,] 9.999994e-01 1.245704e-06 6.228518e-07
[20,] 1.000000e+00 3.410970e-08 1.705485e-08
[21,] 1.000000e+00 1.769866e-08 8.849328e-09
[22,] 1.000000e+00 3.520721e-08 1.760360e-08
[23,] 1.000000e+00 8.692367e-08 4.346184e-08
[24,] 9.999999e-01 2.294003e-07 1.147001e-07
[25,] 9.999998e-01 4.880452e-07 2.440226e-07
[26,] 9.999994e-01 1.260031e-06 6.300156e-07
[27,] 9.999988e-01 2.380729e-06 1.190364e-06
[28,] 9.999989e-01 2.219691e-06 1.109846e-06
[29,] 9.999998e-01 3.716724e-07 1.858362e-07
[30,] 1.000000e+00 5.646232e-12 2.823116e-12
[31,] 1.000000e+00 5.213006e-12 2.606503e-12
[32,] 1.000000e+00 1.641768e-11 8.208838e-12
[33,] 1.000000e+00 5.538983e-11 2.769492e-11
[34,] 1.000000e+00 1.376863e-10 6.884315e-11
[35,] 1.000000e+00 3.730016e-10 1.865008e-10
[36,] 1.000000e+00 7.997489e-10 3.998744e-10
[37,] 1.000000e+00 2.960320e-09 1.480160e-09
[38,] 1.000000e+00 9.494406e-09 4.747203e-09
[39,] 1.000000e+00 8.456122e-09 4.228061e-09
[40,] 1.000000e+00 1.795674e-09 8.978369e-10
[41,] 1.000000e+00 7.116255e-11 3.558127e-11
[42,] 1.000000e+00 3.139012e-14 1.569506e-14
[43,] 1.000000e+00 5.758598e-15 2.879299e-15
[44,] 1.000000e+00 3.114019e-14 1.557009e-14
[45,] 1.000000e+00 2.653077e-13 1.326538e-13
[46,] 1.000000e+00 1.925480e-12 9.627400e-13
[47,] 1.000000e+00 1.334327e-11 6.671634e-12
[48,] 1.000000e+00 9.848228e-11 4.924114e-11
[49,] 1.000000e+00 7.866014e-11 3.933007e-11
[50,] 1.000000e+00 1.727153e-10 8.635766e-11
[51,] 1.000000e+00 2.223358e-09 1.111679e-09
[52,] 1.000000e+00 2.708328e-08 1.354164e-08
[53,] 9.999998e-01 3.040531e-07 1.520265e-07
[54,] 9.999996e-01 7.233939e-07 3.616970e-07
[55,] 9.999995e-01 1.068106e-06 5.340531e-07
[56,] 9.999965e-01 6.921798e-06 3.460899e-06
[57,] 9.999153e-01 1.694371e-04 8.471854e-05
> postscript(file="/var/fisher/rcomp/tmp/1yhif1352144493.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> plot(x[,1], type='l', main='Actuals and Interpolation', ylab='value of Actuals and Interpolation (dots)', xlab='time or index')
> points(x[,1]-mysum$resid)
> grid()
> dev.off()
null device
1
> postscript(file="/var/fisher/rcomp/tmp/2ndrv1352144493.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> plot(mysum$resid, type='b', pch=19, main='Residuals', ylab='value of Residuals', xlab='time or index')
> grid()
> dev.off()
null device
1
> postscript(file="/var/fisher/rcomp/tmp/3vufe1352144493.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> hist(mysum$resid, main='Residual Histogram', xlab='values of Residuals')
> grid()
> dev.off()
null device
1
> postscript(file="/var/fisher/rcomp/tmp/4d3dl1352144493.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> densityplot(~mysum$resid,col='black',main='Residual Density Plot', xlab='values of Residuals')
> dev.off()
null device
1
> postscript(file="/var/fisher/rcomp/tmp/5n0bb1352144493.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> qqnorm(mysum$resid, main='Residual Normal Q-Q Plot')
> qqline(mysum$resid)
> grid()
> dev.off()
null device
1
> (myerror <- as.ts(mysum$resid))
Time Series:
Start = 1
End = 70
Frequency = 1
1 2 3 4 5 6 7
-227.8248 -1049.0208 -1565.0344 -1706.8722 -1556.4871 -1279.3915 303.2590
8 9 10 11 12 13 14
-1109.8581 -2047.0359 -716.2570 -1770.6364 -2023.9906 1348.1733 1592.4001
15 16 17 18 19 20 21
661.9249 -447.7137 -771.4601 -1369.6091 -343.7758 -1218.6103 -2668.4878
22 23 24 25 26 27 28
-2490.1296 -2903.0178 -3387.3351 3957.7392 4398.1940 4054.7597 3350.1194
29 30 31 32 33 34 35
2510.1045 1826.0872 3181.1306 2092.3711 1006.2658 755.4728 -367.2035
36 37 38 39 40 41 42
-1493.6332 1621.9003 4456.5933 2752.3898 1506.5362 854.5266 173.1119
43 44 45 46 47 48 49
1687.9648 1722.9625 1988.1997 1931.9543 1267.1139 568.8735 -910.3478
50 51 52 53 54 55 56
-223.3749 -704.5735 -950.8032 -1129.5335 -1335.8260 -922.0522 -1373.5961
57 58 59 60 61 62 63
-2213.4009 -1999.6508 -2280.4024 -2325.6691 141.3341 688.2326 620.4805
64 65 66 67 68 69 70
130.4141 -786.0677 -855.6553 -506.7264 -389.9802 -916.6877 -812.8575
> postscript(file="/var/fisher/rcomp/tmp/6qlqm1352144493.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> dum <- cbind(lag(myerror,k=1),myerror)
> dum
Time Series:
Start = 0
End = 70
Frequency = 1
lag(myerror, k = 1) myerror
0 -227.8248 NA
1 -1049.0208 -227.8248
2 -1565.0344 -1049.0208
3 -1706.8722 -1565.0344
4 -1556.4871 -1706.8722
5 -1279.3915 -1556.4871
6 303.2590 -1279.3915
7 -1109.8581 303.2590
8 -2047.0359 -1109.8581
9 -716.2570 -2047.0359
10 -1770.6364 -716.2570
11 -2023.9906 -1770.6364
12 1348.1733 -2023.9906
13 1592.4001 1348.1733
14 661.9249 1592.4001
15 -447.7137 661.9249
16 -771.4601 -447.7137
17 -1369.6091 -771.4601
18 -343.7758 -1369.6091
19 -1218.6103 -343.7758
20 -2668.4878 -1218.6103
21 -2490.1296 -2668.4878
22 -2903.0178 -2490.1296
23 -3387.3351 -2903.0178
24 3957.7392 -3387.3351
25 4398.1940 3957.7392
26 4054.7597 4398.1940
27 3350.1194 4054.7597
28 2510.1045 3350.1194
29 1826.0872 2510.1045
30 3181.1306 1826.0872
31 2092.3711 3181.1306
32 1006.2658 2092.3711
33 755.4728 1006.2658
34 -367.2035 755.4728
35 -1493.6332 -367.2035
36 1621.9003 -1493.6332
37 4456.5933 1621.9003
38 2752.3898 4456.5933
39 1506.5362 2752.3898
40 854.5266 1506.5362
41 173.1119 854.5266
42 1687.9648 173.1119
43 1722.9625 1687.9648
44 1988.1997 1722.9625
45 1931.9543 1988.1997
46 1267.1139 1931.9543
47 568.8735 1267.1139
48 -910.3478 568.8735
49 -223.3749 -910.3478
50 -704.5735 -223.3749
51 -950.8032 -704.5735
52 -1129.5335 -950.8032
53 -1335.8260 -1129.5335
54 -922.0522 -1335.8260
55 -1373.5961 -922.0522
56 -2213.4009 -1373.5961
57 -1999.6508 -2213.4009
58 -2280.4024 -1999.6508
59 -2325.6691 -2280.4024
60 141.3341 -2325.6691
61 688.2326 141.3341
62 620.4805 688.2326
63 130.4141 620.4805
64 -786.0677 130.4141
65 -855.6553 -786.0677
66 -506.7264 -855.6553
67 -389.9802 -506.7264
68 -916.6877 -389.9802
69 -812.8575 -916.6877
70 NA -812.8575
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -1049.0208 -227.8248
[2,] -1565.0344 -1049.0208
[3,] -1706.8722 -1565.0344
[4,] -1556.4871 -1706.8722
[5,] -1279.3915 -1556.4871
[6,] 303.2590 -1279.3915
[7,] -1109.8581 303.2590
[8,] -2047.0359 -1109.8581
[9,] -716.2570 -2047.0359
[10,] -1770.6364 -716.2570
[11,] -2023.9906 -1770.6364
[12,] 1348.1733 -2023.9906
[13,] 1592.4001 1348.1733
[14,] 661.9249 1592.4001
[15,] -447.7137 661.9249
[16,] -771.4601 -447.7137
[17,] -1369.6091 -771.4601
[18,] -343.7758 -1369.6091
[19,] -1218.6103 -343.7758
[20,] -2668.4878 -1218.6103
[21,] -2490.1296 -2668.4878
[22,] -2903.0178 -2490.1296
[23,] -3387.3351 -2903.0178
[24,] 3957.7392 -3387.3351
[25,] 4398.1940 3957.7392
[26,] 4054.7597 4398.1940
[27,] 3350.1194 4054.7597
[28,] 2510.1045 3350.1194
[29,] 1826.0872 2510.1045
[30,] 3181.1306 1826.0872
[31,] 2092.3711 3181.1306
[32,] 1006.2658 2092.3711
[33,] 755.4728 1006.2658
[34,] -367.2035 755.4728
[35,] -1493.6332 -367.2035
[36,] 1621.9003 -1493.6332
[37,] 4456.5933 1621.9003
[38,] 2752.3898 4456.5933
[39,] 1506.5362 2752.3898
[40,] 854.5266 1506.5362
[41,] 173.1119 854.5266
[42,] 1687.9648 173.1119
[43,] 1722.9625 1687.9648
[44,] 1988.1997 1722.9625
[45,] 1931.9543 1988.1997
[46,] 1267.1139 1931.9543
[47,] 568.8735 1267.1139
[48,] -910.3478 568.8735
[49,] -223.3749 -910.3478
[50,] -704.5735 -223.3749
[51,] -950.8032 -704.5735
[52,] -1129.5335 -950.8032
[53,] -1335.8260 -1129.5335
[54,] -922.0522 -1335.8260
[55,] -1373.5961 -922.0522
[56,] -2213.4009 -1373.5961
[57,] -1999.6508 -2213.4009
[58,] -2280.4024 -1999.6508
[59,] -2325.6691 -2280.4024
[60,] 141.3341 -2325.6691
[61,] 688.2326 141.3341
[62,] 620.4805 688.2326
[63,] 130.4141 620.4805
[64,] -786.0677 130.4141
[65,] -855.6553 -786.0677
[66,] -506.7264 -855.6553
[67,] -389.9802 -506.7264
[68,] -916.6877 -389.9802
[69,] -812.8575 -916.6877
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -1049.0208 -227.8248
2 -1565.0344 -1049.0208
3 -1706.8722 -1565.0344
4 -1556.4871 -1706.8722
5 -1279.3915 -1556.4871
6 303.2590 -1279.3915
7 -1109.8581 303.2590
8 -2047.0359 -1109.8581
9 -716.2570 -2047.0359
10 -1770.6364 -716.2570
11 -2023.9906 -1770.6364
12 1348.1733 -2023.9906
13 1592.4001 1348.1733
14 661.9249 1592.4001
15 -447.7137 661.9249
16 -771.4601 -447.7137
17 -1369.6091 -771.4601
18 -343.7758 -1369.6091
19 -1218.6103 -343.7758
20 -2668.4878 -1218.6103
21 -2490.1296 -2668.4878
22 -2903.0178 -2490.1296
23 -3387.3351 -2903.0178
24 3957.7392 -3387.3351
25 4398.1940 3957.7392
26 4054.7597 4398.1940
27 3350.1194 4054.7597
28 2510.1045 3350.1194
29 1826.0872 2510.1045
30 3181.1306 1826.0872
31 2092.3711 3181.1306
32 1006.2658 2092.3711
33 755.4728 1006.2658
34 -367.2035 755.4728
35 -1493.6332 -367.2035
36 1621.9003 -1493.6332
37 4456.5933 1621.9003
38 2752.3898 4456.5933
39 1506.5362 2752.3898
40 854.5266 1506.5362
41 173.1119 854.5266
42 1687.9648 173.1119
43 1722.9625 1687.9648
44 1988.1997 1722.9625
45 1931.9543 1988.1997
46 1267.1139 1931.9543
47 568.8735 1267.1139
48 -910.3478 568.8735
49 -223.3749 -910.3478
50 -704.5735 -223.3749
51 -950.8032 -704.5735
52 -1129.5335 -950.8032
53 -1335.8260 -1129.5335
54 -922.0522 -1335.8260
55 -1373.5961 -922.0522
56 -2213.4009 -1373.5961
57 -1999.6508 -2213.4009
58 -2280.4024 -1999.6508
59 -2325.6691 -2280.4024
60 141.3341 -2325.6691
61 688.2326 141.3341
62 620.4805 688.2326
63 130.4141 620.4805
64 -786.0677 130.4141
65 -855.6553 -786.0677
66 -506.7264 -855.6553
67 -389.9802 -506.7264
68 -916.6877 -389.9802
69 -812.8575 -916.6877
> plot(z,main=paste('Residual Lag plot, lowess, and regression line'), ylab='values of Residuals', xlab='lagged values of Residuals')
> lines(lowess(z))
> abline(lm(z))
> grid()
> dev.off()
null device
1
> postscript(file="/var/fisher/rcomp/tmp/7yz9o1352144493.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> acf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Autocorrelation Function')
> grid()
> dev.off()
null device
1
> postscript(file="/var/fisher/rcomp/tmp/8j3ag1352144493.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> pacf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Partial Autocorrelation Function')
> grid()
> dev.off()
null device
1
> postscript(file="/var/fisher/rcomp/tmp/945p91352144493.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> opar <- par(mfrow = c(2,2), oma = c(0, 0, 1.1, 0))
> plot(mylm, las = 1, sub='Residual Diagnostics')
> par(opar)
> dev.off()
null device
1
> if (n > n25) {
+ postscript(file="/var/fisher/rcomp/tmp/10kv1f1352144493.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
+ plot(kp3:nmkm3,gqarr[,2], main='Goldfeld-Quandt test',ylab='2-sided p-value',xlab='breakpoint')
+ grid()
+ dev.off()
+ }
null device
1
>
> #Note: the /var/fisher/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/fisher/rcomp/createtable")
>
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Estimated Regression Equation', 1, TRUE)
> a<-table.row.end(a)
> myeq <- colnames(x)[1]
> myeq <- paste(myeq, '[t] = ', sep='')
> for (i in 1:k){
+ if (mysum$coefficients[i,1] > 0) myeq <- paste(myeq, '+', '')
+ myeq <- paste(myeq, mysum$coefficients[i,1], sep=' ')
+ if (rownames(mysum$coefficients)[i] != '(Intercept)') {
+ myeq <- paste(myeq, rownames(mysum$coefficients)[i], sep='')
+ if (rownames(mysum$coefficients)[i] != 't') myeq <- paste(myeq, '[t]', sep='')
+ }
+ }
> myeq <- paste(myeq, ' + e[t]')
> a<-table.row.start(a)
> a<-table.element(a, myeq)
> a<-table.row.end(a)
> a<-table.end(a)
> table.save(a,file="/var/fisher/rcomp/tmp/11rab11352144493.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a,hyperlink('http://www.xycoon.com/ols1.htm','Multiple Linear Regression - Ordinary Least Squares',''), 6, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a,'Variable',header=TRUE)
> a<-table.element(a,'Parameter',header=TRUE)
> a<-table.element(a,'S.D.',header=TRUE)
> a<-table.element(a,'T-STAT
H0: parameter = 0',header=TRUE)
> a<-table.element(a,'2-tail p-value',header=TRUE)
> a<-table.element(a,'1-tail p-value',header=TRUE)
> a<-table.row.end(a)
> for (i in 1:k){
+ a<-table.row.start(a)
+ a<-table.element(a,rownames(mysum$coefficients)[i],header=TRUE)
+ a<-table.element(a,mysum$coefficients[i,1])
+ a<-table.element(a, round(mysum$coefficients[i,2],6))
+ a<-table.element(a, round(mysum$coefficients[i,3],4))
+ a<-table.element(a, round(mysum$coefficients[i,4],6))
+ a<-table.element(a, round(mysum$coefficients[i,4]/2,6))
+ a<-table.row.end(a)
+ }
> a<-table.end(a)
> table.save(a,file="/var/fisher/rcomp/tmp/12x8bk1352144493.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Regression Statistics', 2, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple R',1,TRUE)
> a<-table.element(a, sqrt(mysum$r.squared))
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'R-squared',1,TRUE)
> a<-table.element(a, mysum$r.squared)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Adjusted R-squared',1,TRUE)
> a<-table.element(a, mysum$adj.r.squared)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (value)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[1])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (DF numerator)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[2])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (DF denominator)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[3])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'p-value',1,TRUE)
> a<-table.element(a, 1-pf(mysum$fstatistic[1],mysum$fstatistic[2],mysum$fstatistic[3]))
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Residual Statistics', 2, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Residual Standard Deviation',1,TRUE)
> a<-table.element(a, mysum$sigma)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Sum Squared Residuals',1,TRUE)
> a<-table.element(a, sum(myerror*myerror))
> a<-table.row.end(a)
> a<-table.end(a)
> table.save(a,file="/var/fisher/rcomp/tmp/13cptk1352144493.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Actuals, Interpolation, and Residuals', 4, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Time or Index', 1, TRUE)
> a<-table.element(a, 'Actuals', 1, TRUE)
> a<-table.element(a, 'Interpolation
Forecast', 1, TRUE)
> a<-table.element(a, 'Residuals
Prediction Error', 1, TRUE)
> a<-table.row.end(a)
> for (i in 1:n) {
+ a<-table.row.start(a)
+ a<-table.element(a,i, 1, TRUE)
+ a<-table.element(a,x[i])
+ a<-table.element(a,x[i]-mysum$resid[i])
+ a<-table.element(a,mysum$resid[i])
+ a<-table.row.end(a)
+ }
> a<-table.end(a)
> table.save(a,file="/var/fisher/rcomp/tmp/149vng1352144493.tab")
> if (n > n25) {
+ a<-table.start()
+ a<-table.row.start(a)
+ a<-table.element(a,'Goldfeld-Quandt test for Heteroskedasticity',4,TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'p-values',header=TRUE)
+ a<-table.element(a,'Alternative Hypothesis',3,header=TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'breakpoint index',header=TRUE)
+ a<-table.element(a,'greater',header=TRUE)
+ a<-table.element(a,'2-sided',header=TRUE)
+ a<-table.element(a,'less',header=TRUE)
+ a<-table.row.end(a)
+ for (mypoint in kp3:nmkm3) {
+ a<-table.row.start(a)
+ a<-table.element(a,mypoint,header=TRUE)
+ a<-table.element(a,gqarr[mypoint-kp3+1,1])
+ a<-table.element(a,gqarr[mypoint-kp3+1,2])
+ a<-table.element(a,gqarr[mypoint-kp3+1,3])
+ a<-table.row.end(a)
+ }
+ a<-table.end(a)
+ table.save(a,file="/var/fisher/rcomp/tmp/15kww41352144493.tab")
+ a<-table.start()
+ a<-table.row.start(a)
+ a<-table.element(a,'Meta Analysis of Goldfeld-Quandt test for Heteroskedasticity',4,TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'Description',header=TRUE)
+ a<-table.element(a,'# significant tests',header=TRUE)
+ a<-table.element(a,'% significant tests',header=TRUE)
+ a<-table.element(a,'OK/NOK',header=TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'1% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant1)
+ a<-table.element(a,numsignificant1/numgqtests)
+ if (numsignificant1/numgqtests < 0.01) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'5% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant5)
+ a<-table.element(a,numsignificant5/numgqtests)
+ if (numsignificant5/numgqtests < 0.05) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'10% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant10)
+ a<-table.element(a,numsignificant10/numgqtests)
+ if (numsignificant10/numgqtests < 0.1) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.end(a)
+ table.save(a,file="/var/fisher/rcomp/tmp/16dzzg1352144493.tab")
+ }
>
> try(system("convert tmp/1yhif1352144493.ps tmp/1yhif1352144493.png",intern=TRUE))
character(0)
> try(system("convert tmp/2ndrv1352144493.ps tmp/2ndrv1352144493.png",intern=TRUE))
character(0)
> try(system("convert tmp/3vufe1352144493.ps tmp/3vufe1352144493.png",intern=TRUE))
character(0)
> try(system("convert tmp/4d3dl1352144493.ps tmp/4d3dl1352144493.png",intern=TRUE))
character(0)
> try(system("convert tmp/5n0bb1352144493.ps tmp/5n0bb1352144493.png",intern=TRUE))
character(0)
> try(system("convert tmp/6qlqm1352144493.ps tmp/6qlqm1352144493.png",intern=TRUE))
character(0)
> try(system("convert tmp/7yz9o1352144493.ps tmp/7yz9o1352144493.png",intern=TRUE))
character(0)
> try(system("convert tmp/8j3ag1352144493.ps tmp/8j3ag1352144493.png",intern=TRUE))
character(0)
> try(system("convert tmp/945p91352144493.ps tmp/945p91352144493.png",intern=TRUE))
character(0)
> try(system("convert tmp/10kv1f1352144493.ps tmp/10kv1f1352144493.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
6.092 1.114 7.210