R version 2.9.0 (2009-04-17)
Copyright (C) 2009 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
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(1.579
+ ,9.769
+ ,2.146
+ ,9.321
+ ,2.462
+ ,9.939
+ ,3.695
+ ,9.336
+ ,4.831
+ ,10.195
+ ,5.134
+ ,9.464
+ ,6.250
+ ,10.010
+ ,5.760
+ ,10.213
+ ,6.249
+ ,9.563
+ ,2.917
+ ,9.890
+ ,1.741
+ ,9.305
+ ,2.359
+ ,9.391
+ ,1.511
+ ,9.928
+ ,2.059
+ ,8.686
+ ,2.635
+ ,9.843
+ ,2.867
+ ,9.627
+ ,4.403
+ ,10.074
+ ,5.720
+ ,9.503
+ ,4.502
+ ,10.119
+ ,5.749
+ ,10.000
+ ,5.627
+ ,9.313
+ ,2.846
+ ,9.866
+ ,1.762
+ ,9.172
+ ,2.429
+ ,9.241
+ ,1.169
+ ,9.659
+ ,2.154
+ ,8.904
+ ,2.249
+ ,9.755
+ ,2.687
+ ,9.080
+ ,4.359
+ ,9.435
+ ,5.382
+ ,8.971
+ ,4.459
+ ,10.063
+ ,6.398
+ ,9.793
+ ,4.596
+ ,9.454
+ ,3.024
+ ,9.759
+ ,1.887
+ ,8.820
+ ,2.070
+ ,9.403
+ ,1.351
+ ,9.676
+ ,2.218
+ ,8.642
+ ,2.461
+ ,9.402
+ ,3.028
+ ,9.610
+ ,4.784
+ ,9.294
+ ,4.975
+ ,9.448
+ ,4.607
+ ,10.319
+ ,6.249
+ ,9.548
+ ,4.809
+ ,9.801
+ ,3.157
+ ,9.596
+ ,1.910
+ ,8.923
+ ,2.228
+ ,9.746
+ ,1.594
+ ,9.829
+ ,2.467
+ ,9.125
+ ,2.222
+ ,9.782
+ ,3.607
+ ,9.441
+ ,4.685
+ ,9.162
+ ,4.962
+ ,9.915
+ ,5.770
+ ,10.444
+ ,5.480
+ ,10.209
+ ,5.000
+ ,9.985
+ ,3.228
+ ,9.842
+ ,1.993
+ ,9.429
+ ,2.288
+ ,10.132
+ ,1.580
+ ,9.849
+ ,2.111
+ ,9.172
+ ,2.192
+ ,10.313
+ ,3.601
+ ,9.819
+ ,4.665
+ ,9.955
+ ,4.876
+ ,10.048
+ ,5.813
+ ,10.082
+ ,5.589
+ ,10.541
+ ,5.331
+ ,10.208
+ ,3.075
+ ,10.233
+ ,2.002
+ ,9.439
+ ,2.306
+ ,9.963
+ ,1.507
+ ,10.158
+ ,1.992
+ ,9.225
+ ,2.487
+ ,10.474
+ ,3.490
+ ,9.757
+ ,4.647
+ ,10.490
+ ,5.594
+ ,10.281
+ ,5.611
+ ,10.444
+ ,5.788
+ ,10.640
+ ,6.204
+ ,10.695
+ ,3.013
+ ,10.786
+ ,1.931
+ ,9.832
+ ,2.549
+ ,9.747
+ ,1.504
+ ,10.411
+ ,2.090
+ ,9.511
+ ,2.702
+ ,10.402
+ ,2.939
+ ,9.701
+ ,4.500
+ ,10.540
+ ,6.208
+ ,10.112
+ ,6.415
+ ,10.915
+ ,5.657
+ ,11.183
+ ,5.964
+ ,10.384
+ ,3.163
+ ,10.834
+ ,1.997
+ ,9.886
+ ,2.422
+ ,10.216)
+ ,dim=c(2
+ ,96)
+ ,dimnames=list(c('huwelijken'
+ ,'geboortes')
+ ,1:96))
> y <- array(NA,dim=c(2,96),dimnames=list(c('huwelijken','geboortes'),1:96))
> 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'
> #'GNU S' R Code compiled by R2WASP v. 1.0.44 ()
> #Author: Prof. Dr. P. Wessa
> #To cite this work: AUTHOR(S), (YEAR), YOUR SOFTWARE TITLE (vNUMBER) in Free Statistics Software (v$_version), Office for Research Development and Education, URL http://www.wessa.net/rwasp_YOURPAGE.wasp/
> #Source of accompanying publication: Office for Research, Development, and Education
> #Technical description: Write here your technical program description (don't use hard returns!)
> library(lattice)
> library(lmtest)
Loading required package: zoo
Attaching package: 'zoo'
The following object(s) are masked from package:base :
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
huwelijken geboortes
1 1.579 9.769
2 2.146 9.321
3 2.462 9.939
4 3.695 9.336
5 4.831 10.195
6 5.134 9.464
7 6.250 10.010
8 5.760 10.213
9 6.249 9.563
10 2.917 9.890
11 1.741 9.305
12 2.359 9.391
13 1.511 9.928
14 2.059 8.686
15 2.635 9.843
16 2.867 9.627
17 4.403 10.074
18 5.720 9.503
19 4.502 10.119
20 5.749 10.000
21 5.627 9.313
22 2.846 9.866
23 1.762 9.172
24 2.429 9.241
25 1.169 9.659
26 2.154 8.904
27 2.249 9.755
28 2.687 9.080
29 4.359 9.435
30 5.382 8.971
31 4.459 10.063
32 6.398 9.793
33 4.596 9.454
34 3.024 9.759
35 1.887 8.820
36 2.070 9.403
37 1.351 9.676
38 2.218 8.642
39 2.461 9.402
40 3.028 9.610
41 4.784 9.294
42 4.975 9.448
43 4.607 10.319
44 6.249 9.548
45 4.809 9.801
46 3.157 9.596
47 1.910 8.923
48 2.228 9.746
49 1.594 9.829
50 2.467 9.125
51 2.222 9.782
52 3.607 9.441
53 4.685 9.162
54 4.962 9.915
55 5.770 10.444
56 5.480 10.209
57 5.000 9.985
58 3.228 9.842
59 1.993 9.429
60 2.288 10.132
61 1.580 9.849
62 2.111 9.172
63 2.192 10.313
64 3.601 9.819
65 4.665 9.955
66 4.876 10.048
67 5.813 10.082
68 5.589 10.541
69 5.331 10.208
70 3.075 10.233
71 2.002 9.439
72 2.306 9.963
73 1.507 10.158
74 1.992 9.225
75 2.487 10.474
76 3.490 9.757
77 4.647 10.490
78 5.594 10.281
79 5.611 10.444
80 5.788 10.640
81 6.204 10.695
82 3.013 10.786
83 1.931 9.832
84 2.549 9.747
85 1.504 10.411
86 2.090 9.511
87 2.702 10.402
88 2.939 9.701
89 4.500 10.540
90 6.208 10.112
91 6.415 10.915
92 5.657 11.183
93 5.964 10.384
94 3.163 10.834
95 1.997 9.886
96 2.422 10.216
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) geboortes
-8.060 1.188
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-2.8091 -1.1540 -0.2608 1.2435 2.9616
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -8.060 2.919 -2.761 0.006930 **
geboortes 1.188 0.297 4.002 0.000125 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.485 on 94 degrees of freedom
Multiple R-squared: 0.1456, Adjusted R-squared: 0.1365
F-statistic: 16.01 on 1 and 94 DF, p-value: 0.0001253
> 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,] 0.5420045 0.9159909 0.4579955
[2,] 0.6845747 0.6308506 0.3154253
[3,] 0.8034983 0.3930035 0.1965017
[4,] 0.7519272 0.4961456 0.2480728
[5,] 0.8611520 0.2776960 0.1388480
[6,] 0.8420256 0.3159488 0.1579744
[7,] 0.8281321 0.3437358 0.1718679
[8,] 0.7769357 0.4461285 0.2230643
[9,] 0.8630958 0.2738085 0.1369042
[10,] 0.8108315 0.3783371 0.1891685
[11,] 0.7793056 0.4413888 0.2206944
[12,] 0.7201779 0.5596442 0.2798221
[13,] 0.6524483 0.6951033 0.3475517
[14,] 0.7554840 0.4890321 0.2445160
[15,] 0.6948875 0.6102249 0.3051125
[16,] 0.7070149 0.5859701 0.2929851
[17,] 0.7988036 0.4023927 0.2011964
[18,] 0.7717735 0.4564531 0.2282265
[19,] 0.7491011 0.5017978 0.2508989
[20,] 0.6982392 0.6035216 0.3017608
[21,] 0.7726069 0.4547863 0.2273931
[22,] 0.7202240 0.5595519 0.2797760
[23,] 0.7082186 0.5835628 0.2917814
[24,] 0.6494330 0.7011340 0.3505670
[25,] 0.6261730 0.7476541 0.3738270
[26,] 0.7553972 0.4892056 0.2446028
[27,] 0.7083661 0.5832679 0.2916339
[28,] 0.8142508 0.3714984 0.1857492
[29,] 0.8056589 0.3886821 0.1943411
[30,] 0.7692660 0.4614679 0.2307340
[31,] 0.7277210 0.5445580 0.2722790
[32,] 0.7034083 0.5931833 0.2965917
[33,] 0.7546538 0.4906924 0.2453462
[34,] 0.7050801 0.5898397 0.2949199
[35,] 0.6617165 0.6765671 0.3382835
[36,] 0.6085735 0.7828530 0.3914265
[37,] 0.6357325 0.7285350 0.3642675
[38,] 0.6644646 0.6710709 0.3355354
[39,] 0.6114428 0.7771143 0.3885572
[40,] 0.7736001 0.4527998 0.2263999
[41,] 0.7628061 0.4743879 0.2371939
[42,] 0.7178313 0.5643373 0.2821687
[43,] 0.6737901 0.6524199 0.3262099
[44,] 0.6591341 0.6817319 0.3408659
[45,] 0.7005949 0.5988102 0.2994051
[46,] 0.6503438 0.6993124 0.3496562
[47,] 0.6351601 0.7296798 0.3648399
[48,] 0.5904884 0.8190231 0.4095116
[49,] 0.6665111 0.6669778 0.3334889
[50,] 0.6644038 0.6711923 0.3355962
[51,] 0.6566148 0.6867703 0.3433852
[52,] 0.6579430 0.6841141 0.3420570
[53,] 0.6555350 0.6889299 0.3444650
[54,] 0.6037955 0.7924090 0.3962045
[55,] 0.5641299 0.8717401 0.4358701
[56,] 0.5755198 0.8489604 0.4244802
[57,] 0.6077247 0.7845506 0.3922753
[58,] 0.5543218 0.8913564 0.4456782
[59,] 0.5991507 0.8016987 0.4008493
[60,] 0.5415624 0.9168752 0.4584376
[61,] 0.5188929 0.9622143 0.4811071
[62,] 0.5021033 0.9957934 0.4978967
[63,] 0.5810246 0.8379509 0.4189754
[64,] 0.5564212 0.8871575 0.4435788
[65,] 0.5667674 0.8664652 0.4332326
[66,] 0.5204745 0.9590511 0.4795255
[67,] 0.4653758 0.9307516 0.5346242
[68,] 0.4333777 0.8667554 0.5666223
[69,] 0.5188186 0.9623628 0.4811814
[70,] 0.4551063 0.9102125 0.5448937
[71,] 0.4907664 0.9815329 0.5092336
[72,] 0.4287776 0.8575551 0.5712224
[73,] 0.3587474 0.7174948 0.6412526
[74,] 0.3762651 0.7525302 0.6237349
[75,] 0.3704927 0.7409854 0.6295073
[76,] 0.3523524 0.7047048 0.6476476
[77,] 0.3806257 0.7612514 0.6193743
[78,] 0.3898895 0.7797791 0.6101105
[79,] 0.3417068 0.6834135 0.6582932
[80,] 0.2657257 0.5314514 0.7342743
[81,] 0.4203608 0.8407216 0.5796392
[82,] 0.3321482 0.6642964 0.6678518
[83,] 0.3279222 0.6558444 0.6720778
[84,] 0.2322344 0.4644689 0.7677656
[85,] 0.1495994 0.2991987 0.8504006
[86,] 0.3116615 0.6233230 0.6883385
[87,] 0.2769724 0.5539448 0.7230276
> postscript(file="/var/www/html/rcomp/tmp/1vljv1290938526.ps",horizontal=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/www/html/rcomp/tmp/2odiy1290938526.ps",horizontal=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/www/html/rcomp/tmp/3odiy1290938526.ps",horizontal=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/www/html/rcomp/tmp/4h4z11290938526.ps",horizontal=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/www/html/rcomp/tmp/5h4z11290938526.ps",horizontal=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 = 96
Frequency = 1
1 2 3 4 5 6
-1.971047286 -0.871599326 -1.290092270 0.659573175 0.774651753 1.946445187
7 8 9 10 11 12
2.413524236 1.682258755 2.943783696 -0.776855775 -1.257583328 -0.741794320
13 14 15 16 17 18
-2.228018771 -0.203901884 -1.002996279 -0.514280298 0.490460242 2.486093690
19 20 21 22 23 24
0.535977746 1.924409235 2.618908673 -0.819331777 -1.078512840 -0.493519334
25 26 27 28 29 30
-2.250312296 -0.367994864 -1.284408287 -0.044170848 1.205911684 2.780375642
31 32 33 34 35 36
0.559533741 2.819428717 1.420330186 -0.514162287 -0.535160872 -1.045056319
37 38 39 40 41 42
-2.088516794 0.007392112 -0.652867819 -0.333075800 1.798490171 1.806461185
43 44 45 46 47 48
0.403277764 2.961611194 1.220920717 -0.187436801 -0.634576362 -1.294711788
49 50 51 52 53 54
-2.027357280 -0.317653344 -1.343497784 0.446780685 1.856372159 1.238431728
55 56 57 58 59 60
1.417715276 1.407012754 1.193236734 -0.408807779 -1.152957316 -1.693472753
61 62 63 64 65 66
-2.065127278 -0.729512840 -2.004591236 -0.008472281 0.893891731 0.994361240
67 68 69 70 71 72
1.890952243 1.121430785 1.259201254 -1.026511243 -1.155842316 -1.474616268
73 74 75 76 77 78
-2.505373750 -0.911503335 -1.900939722 -0.045785287 0.240044280 1.435440761
79 80 81 82 83 84
1.258715276 1.202769294 1.553401799 -1.745751693 -1.693922780 -0.974900288
85 86 87 88 89 90
-2.809064227 -1.153414309 -1.600367728 -0.530229292 0.033619284 2.250297246
91 92 93 94 95 96
1.502931819 0.426413843 1.683025270 -1.652799689 -1.692101775 -1.659306745
> postscript(file="/var/www/html/rcomp/tmp/6h4z11290938526.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> dum <- cbind(lag(myerror,k=1),myerror)
> dum
Time Series:
Start = 0
End = 96
Frequency = 1
lag(myerror, k = 1) myerror
0 -1.971047286 NA
1 -0.871599326 -1.971047286
2 -1.290092270 -0.871599326
3 0.659573175 -1.290092270
4 0.774651753 0.659573175
5 1.946445187 0.774651753
6 2.413524236 1.946445187
7 1.682258755 2.413524236
8 2.943783696 1.682258755
9 -0.776855775 2.943783696
10 -1.257583328 -0.776855775
11 -0.741794320 -1.257583328
12 -2.228018771 -0.741794320
13 -0.203901884 -2.228018771
14 -1.002996279 -0.203901884
15 -0.514280298 -1.002996279
16 0.490460242 -0.514280298
17 2.486093690 0.490460242
18 0.535977746 2.486093690
19 1.924409235 0.535977746
20 2.618908673 1.924409235
21 -0.819331777 2.618908673
22 -1.078512840 -0.819331777
23 -0.493519334 -1.078512840
24 -2.250312296 -0.493519334
25 -0.367994864 -2.250312296
26 -1.284408287 -0.367994864
27 -0.044170848 -1.284408287
28 1.205911684 -0.044170848
29 2.780375642 1.205911684
30 0.559533741 2.780375642
31 2.819428717 0.559533741
32 1.420330186 2.819428717
33 -0.514162287 1.420330186
34 -0.535160872 -0.514162287
35 -1.045056319 -0.535160872
36 -2.088516794 -1.045056319
37 0.007392112 -2.088516794
38 -0.652867819 0.007392112
39 -0.333075800 -0.652867819
40 1.798490171 -0.333075800
41 1.806461185 1.798490171
42 0.403277764 1.806461185
43 2.961611194 0.403277764
44 1.220920717 2.961611194
45 -0.187436801 1.220920717
46 -0.634576362 -0.187436801
47 -1.294711788 -0.634576362
48 -2.027357280 -1.294711788
49 -0.317653344 -2.027357280
50 -1.343497784 -0.317653344
51 0.446780685 -1.343497784
52 1.856372159 0.446780685
53 1.238431728 1.856372159
54 1.417715276 1.238431728
55 1.407012754 1.417715276
56 1.193236734 1.407012754
57 -0.408807779 1.193236734
58 -1.152957316 -0.408807779
59 -1.693472753 -1.152957316
60 -2.065127278 -1.693472753
61 -0.729512840 -2.065127278
62 -2.004591236 -0.729512840
63 -0.008472281 -2.004591236
64 0.893891731 -0.008472281
65 0.994361240 0.893891731
66 1.890952243 0.994361240
67 1.121430785 1.890952243
68 1.259201254 1.121430785
69 -1.026511243 1.259201254
70 -1.155842316 -1.026511243
71 -1.474616268 -1.155842316
72 -2.505373750 -1.474616268
73 -0.911503335 -2.505373750
74 -1.900939722 -0.911503335
75 -0.045785287 -1.900939722
76 0.240044280 -0.045785287
77 1.435440761 0.240044280
78 1.258715276 1.435440761
79 1.202769294 1.258715276
80 1.553401799 1.202769294
81 -1.745751693 1.553401799
82 -1.693922780 -1.745751693
83 -0.974900288 -1.693922780
84 -2.809064227 -0.974900288
85 -1.153414309 -2.809064227
86 -1.600367728 -1.153414309
87 -0.530229292 -1.600367728
88 0.033619284 -0.530229292
89 2.250297246 0.033619284
90 1.502931819 2.250297246
91 0.426413843 1.502931819
92 1.683025270 0.426413843
93 -1.652799689 1.683025270
94 -1.692101775 -1.652799689
95 -1.659306745 -1.692101775
96 NA -1.659306745
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -0.871599326 -1.971047286
[2,] -1.290092270 -0.871599326
[3,] 0.659573175 -1.290092270
[4,] 0.774651753 0.659573175
[5,] 1.946445187 0.774651753
[6,] 2.413524236 1.946445187
[7,] 1.682258755 2.413524236
[8,] 2.943783696 1.682258755
[9,] -0.776855775 2.943783696
[10,] -1.257583328 -0.776855775
[11,] -0.741794320 -1.257583328
[12,] -2.228018771 -0.741794320
[13,] -0.203901884 -2.228018771
[14,] -1.002996279 -0.203901884
[15,] -0.514280298 -1.002996279
[16,] 0.490460242 -0.514280298
[17,] 2.486093690 0.490460242
[18,] 0.535977746 2.486093690
[19,] 1.924409235 0.535977746
[20,] 2.618908673 1.924409235
[21,] -0.819331777 2.618908673
[22,] -1.078512840 -0.819331777
[23,] -0.493519334 -1.078512840
[24,] -2.250312296 -0.493519334
[25,] -0.367994864 -2.250312296
[26,] -1.284408287 -0.367994864
[27,] -0.044170848 -1.284408287
[28,] 1.205911684 -0.044170848
[29,] 2.780375642 1.205911684
[30,] 0.559533741 2.780375642
[31,] 2.819428717 0.559533741
[32,] 1.420330186 2.819428717
[33,] -0.514162287 1.420330186
[34,] -0.535160872 -0.514162287
[35,] -1.045056319 -0.535160872
[36,] -2.088516794 -1.045056319
[37,] 0.007392112 -2.088516794
[38,] -0.652867819 0.007392112
[39,] -0.333075800 -0.652867819
[40,] 1.798490171 -0.333075800
[41,] 1.806461185 1.798490171
[42,] 0.403277764 1.806461185
[43,] 2.961611194 0.403277764
[44,] 1.220920717 2.961611194
[45,] -0.187436801 1.220920717
[46,] -0.634576362 -0.187436801
[47,] -1.294711788 -0.634576362
[48,] -2.027357280 -1.294711788
[49,] -0.317653344 -2.027357280
[50,] -1.343497784 -0.317653344
[51,] 0.446780685 -1.343497784
[52,] 1.856372159 0.446780685
[53,] 1.238431728 1.856372159
[54,] 1.417715276 1.238431728
[55,] 1.407012754 1.417715276
[56,] 1.193236734 1.407012754
[57,] -0.408807779 1.193236734
[58,] -1.152957316 -0.408807779
[59,] -1.693472753 -1.152957316
[60,] -2.065127278 -1.693472753
[61,] -0.729512840 -2.065127278
[62,] -2.004591236 -0.729512840
[63,] -0.008472281 -2.004591236
[64,] 0.893891731 -0.008472281
[65,] 0.994361240 0.893891731
[66,] 1.890952243 0.994361240
[67,] 1.121430785 1.890952243
[68,] 1.259201254 1.121430785
[69,] -1.026511243 1.259201254
[70,] -1.155842316 -1.026511243
[71,] -1.474616268 -1.155842316
[72,] -2.505373750 -1.474616268
[73,] -0.911503335 -2.505373750
[74,] -1.900939722 -0.911503335
[75,] -0.045785287 -1.900939722
[76,] 0.240044280 -0.045785287
[77,] 1.435440761 0.240044280
[78,] 1.258715276 1.435440761
[79,] 1.202769294 1.258715276
[80,] 1.553401799 1.202769294
[81,] -1.745751693 1.553401799
[82,] -1.693922780 -1.745751693
[83,] -0.974900288 -1.693922780
[84,] -2.809064227 -0.974900288
[85,] -1.153414309 -2.809064227
[86,] -1.600367728 -1.153414309
[87,] -0.530229292 -1.600367728
[88,] 0.033619284 -0.530229292
[89,] 2.250297246 0.033619284
[90,] 1.502931819 2.250297246
[91,] 0.426413843 1.502931819
[92,] 1.683025270 0.426413843
[93,] -1.652799689 1.683025270
[94,] -1.692101775 -1.652799689
[95,] -1.659306745 -1.692101775
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -0.871599326 -1.971047286
2 -1.290092270 -0.871599326
3 0.659573175 -1.290092270
4 0.774651753 0.659573175
5 1.946445187 0.774651753
6 2.413524236 1.946445187
7 1.682258755 2.413524236
8 2.943783696 1.682258755
9 -0.776855775 2.943783696
10 -1.257583328 -0.776855775
11 -0.741794320 -1.257583328
12 -2.228018771 -0.741794320
13 -0.203901884 -2.228018771
14 -1.002996279 -0.203901884
15 -0.514280298 -1.002996279
16 0.490460242 -0.514280298
17 2.486093690 0.490460242
18 0.535977746 2.486093690
19 1.924409235 0.535977746
20 2.618908673 1.924409235
21 -0.819331777 2.618908673
22 -1.078512840 -0.819331777
23 -0.493519334 -1.078512840
24 -2.250312296 -0.493519334
25 -0.367994864 -2.250312296
26 -1.284408287 -0.367994864
27 -0.044170848 -1.284408287
28 1.205911684 -0.044170848
29 2.780375642 1.205911684
30 0.559533741 2.780375642
31 2.819428717 0.559533741
32 1.420330186 2.819428717
33 -0.514162287 1.420330186
34 -0.535160872 -0.514162287
35 -1.045056319 -0.535160872
36 -2.088516794 -1.045056319
37 0.007392112 -2.088516794
38 -0.652867819 0.007392112
39 -0.333075800 -0.652867819
40 1.798490171 -0.333075800
41 1.806461185 1.798490171
42 0.403277764 1.806461185
43 2.961611194 0.403277764
44 1.220920717 2.961611194
45 -0.187436801 1.220920717
46 -0.634576362 -0.187436801
47 -1.294711788 -0.634576362
48 -2.027357280 -1.294711788
49 -0.317653344 -2.027357280
50 -1.343497784 -0.317653344
51 0.446780685 -1.343497784
52 1.856372159 0.446780685
53 1.238431728 1.856372159
54 1.417715276 1.238431728
55 1.407012754 1.417715276
56 1.193236734 1.407012754
57 -0.408807779 1.193236734
58 -1.152957316 -0.408807779
59 -1.693472753 -1.152957316
60 -2.065127278 -1.693472753
61 -0.729512840 -2.065127278
62 -2.004591236 -0.729512840
63 -0.008472281 -2.004591236
64 0.893891731 -0.008472281
65 0.994361240 0.893891731
66 1.890952243 0.994361240
67 1.121430785 1.890952243
68 1.259201254 1.121430785
69 -1.026511243 1.259201254
70 -1.155842316 -1.026511243
71 -1.474616268 -1.155842316
72 -2.505373750 -1.474616268
73 -0.911503335 -2.505373750
74 -1.900939722 -0.911503335
75 -0.045785287 -1.900939722
76 0.240044280 -0.045785287
77 1.435440761 0.240044280
78 1.258715276 1.435440761
79 1.202769294 1.258715276
80 1.553401799 1.202769294
81 -1.745751693 1.553401799
82 -1.693922780 -1.745751693
83 -0.974900288 -1.693922780
84 -2.809064227 -0.974900288
85 -1.153414309 -2.809064227
86 -1.600367728 -1.153414309
87 -0.530229292 -1.600367728
88 0.033619284 -0.530229292
89 2.250297246 0.033619284
90 1.502931819 2.250297246
91 0.426413843 1.502931819
92 1.683025270 0.426413843
93 -1.652799689 1.683025270
94 -1.692101775 -1.652799689
95 -1.659306745 -1.692101775
> 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/www/html/rcomp/tmp/79vgm1290938526.ps",horizontal=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/www/html/rcomp/tmp/89vgm1290938526.ps",horizontal=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/www/html/rcomp/tmp/92mg71290938526.ps",horizontal=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/www/html/rcomp/tmp/102mg71290938526.ps",horizontal=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/www/html/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/html/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/www/html/rcomp/tmp/11o5ed1290938526.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/www/html/rcomp/tmp/12r5d11290938526.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/www/html/rcomp/tmp/13ypad1290938526.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/www/html/rcomp/tmp/14qgrf1290938526.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/www/html/rcomp/tmp/15ug7l1290938526.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/www/html/rcomp/tmp/16q8nu1290938526.tab")
+ }
>
> try(system("convert tmp/1vljv1290938526.ps tmp/1vljv1290938526.png",intern=TRUE))
character(0)
> try(system("convert tmp/2odiy1290938526.ps tmp/2odiy1290938526.png",intern=TRUE))
character(0)
> try(system("convert tmp/3odiy1290938526.ps tmp/3odiy1290938526.png",intern=TRUE))
character(0)
> try(system("convert tmp/4h4z11290938526.ps tmp/4h4z11290938526.png",intern=TRUE))
character(0)
> try(system("convert tmp/5h4z11290938526.ps tmp/5h4z11290938526.png",intern=TRUE))
character(0)
> try(system("convert tmp/6h4z11290938526.ps tmp/6h4z11290938526.png",intern=TRUE))
character(0)
> try(system("convert tmp/79vgm1290938526.ps tmp/79vgm1290938526.png",intern=TRUE))
character(0)
> try(system("convert tmp/89vgm1290938526.ps tmp/89vgm1290938526.png",intern=TRUE))
character(0)
> try(system("convert tmp/92mg71290938526.ps tmp/92mg71290938526.png",intern=TRUE))
character(0)
> try(system("convert tmp/102mg71290938526.ps tmp/102mg71290938526.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
2.872 1.628 6.504