R version 2.8.0 (2008-10-20)
Copyright (C) 2008 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.
Natural language support but running in an English locale
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(9769
+ ,1579
+ ,9321
+ ,2146
+ ,9939
+ ,2462
+ ,9336
+ ,3695
+ ,10195
+ ,4831
+ ,9464
+ ,5134
+ ,10010
+ ,6250
+ ,10213
+ ,5760
+ ,9563
+ ,6249
+ ,9890
+ ,2917
+ ,9305
+ ,1741
+ ,9391
+ ,2359
+ ,9928
+ ,1511
+ ,8686
+ ,2059
+ ,9843
+ ,2635
+ ,9627
+ ,2867
+ ,10074
+ ,4403
+ ,9503
+ ,5720
+ ,10119
+ ,4502
+ ,10000
+ ,5749
+ ,9313
+ ,5627
+ ,9866
+ ,2846
+ ,9172
+ ,1762
+ ,9241
+ ,2429
+ ,9659
+ ,1169
+ ,8904
+ ,2154
+ ,9755
+ ,2249
+ ,9080
+ ,2687
+ ,9435
+ ,4359
+ ,8971
+ ,5382
+ ,10063
+ ,4459
+ ,9793
+ ,6398
+ ,9454
+ ,4596
+ ,9759
+ ,3024
+ ,8820
+ ,1887
+ ,9403
+ ,2070
+ ,9676
+ ,1351
+ ,8642
+ ,2218
+ ,9402
+ ,2461
+ ,9610
+ ,3028
+ ,9294
+ ,4784
+ ,9448
+ ,4975
+ ,10319
+ ,4607
+ ,9548
+ ,6249
+ ,9801
+ ,4809
+ ,9596
+ ,3157
+ ,8923
+ ,1910
+ ,9746
+ ,2228
+ ,9829
+ ,1594
+ ,9125
+ ,2467
+ ,9782
+ ,2222
+ ,9441
+ ,3607
+ ,9162
+ ,4685
+ ,9915
+ ,4962
+ ,10444
+ ,5770
+ ,10209
+ ,5480
+ ,9985
+ ,5000
+ ,9842
+ ,3228
+ ,9429
+ ,1993
+ ,10132
+ ,2288
+ ,9849
+ ,1580
+ ,9172
+ ,2111
+ ,10313
+ ,2192
+ ,9819
+ ,3601
+ ,9955
+ ,4665
+ ,10048
+ ,4876
+ ,10082
+ ,5813
+ ,10541
+ ,5589
+ ,10208
+ ,5331
+ ,10233
+ ,3075
+ ,9439
+ ,2002
+ ,9963
+ ,2306
+ ,10158
+ ,1507
+ ,9225
+ ,1992
+ ,10474
+ ,2487
+ ,9757
+ ,3490
+ ,10490
+ ,4647
+ ,10281
+ ,5594
+ ,10444
+ ,5611
+ ,10640
+ ,5788
+ ,10695
+ ,6204
+ ,10786
+ ,3013
+ ,9832
+ ,1931
+ ,9747
+ ,2549
+ ,10411
+ ,1504
+ ,9511
+ ,2090
+ ,10402
+ ,2702
+ ,9701
+ ,2939
+ ,10540
+ ,4500
+ ,10112
+ ,6208
+ ,10915
+ ,6415
+ ,11183
+ ,5657
+ ,10384
+ ,5964
+ ,10834
+ ,3163
+ ,9886
+ ,1997
+ ,10216
+ ,2422)
+ ,dim=c(2
+ ,96)
+ ,dimnames=list(c('geboortes'
+ ,'huwelijken')
+ ,1:96))
> y <- array(NA,dim=c(2,96),dimnames=list(c('geboortes','huwelijken'),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
geboortes huwelijken
1 9769 1579
2 9321 2146
3 9939 2462
4 9336 3695
5 10195 4831
6 9464 5134
7 10010 6250
8 10213 5760
9 9563 6249
10 9890 2917
11 9305 1741
12 9391 2359
13 9928 1511
14 8686 2059
15 9843 2635
16 9627 2867
17 10074 4403
18 9503 5720
19 10119 4502
20 10000 5749
21 9313 5627
22 9866 2846
23 9172 1762
24 9241 2429
25 9659 1169
26 8904 2154
27 9755 2249
28 9080 2687
29 9435 4359
30 8971 5382
31 10063 4459
32 9793 6398
33 9454 4596
34 9759 3024
35 8820 1887
36 9403 2070
37 9676 1351
38 8642 2218
39 9402 2461
40 9610 3028
41 9294 4784
42 9448 4975
43 10319 4607
44 9548 6249
45 9801 4809
46 9596 3157
47 8923 1910
48 9746 2228
49 9829 1594
50 9125 2467
51 9782 2222
52 9441 3607
53 9162 4685
54 9915 4962
55 10444 5770
56 10209 5480
57 9985 5000
58 9842 3228
59 9429 1993
60 10132 2288
61 9849 1580
62 9172 2111
63 10313 2192
64 9819 3601
65 9955 4665
66 10048 4876
67 10082 5813
68 10541 5589
69 10208 5331
70 10233 3075
71 9439 2002
72 9963 2306
73 10158 1507
74 9225 1992
75 10474 2487
76 9757 3490
77 10490 4647
78 10281 5594
79 10444 5611
80 10640 5788
81 10695 6204
82 10786 3013
83 9832 1931
84 9747 2549
85 10411 1504
86 9511 2090
87 10402 2702
88 9701 2939
89 10540 4500
90 10112 6208
91 10915 6415
92 11183 5657
93 10384 5964
94 10834 3163
95 9886 1997
96 10216 2422
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) huwelijken
9374.9968 0.1225
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-1063.20 -368.18 36.70 269.43 1115.12
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 9.375e+03 1.206e+02 77.715 < 2e-16 ***
huwelijken 1.225e-01 3.061e-02 4.002 0.000125 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 476.6 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.45189814 0.9037963 0.5481019
[2,] 0.38508464 0.7701693 0.6149154
[3,] 0.26870134 0.5374027 0.7312987
[4,] 0.21490342 0.4298068 0.7850966
[5,] 0.19627293 0.3925459 0.8037271
[6,] 0.13664427 0.2732885 0.8633557
[7,] 0.10914452 0.2182890 0.8908555
[8,] 0.07601047 0.1520209 0.9239895
[9,] 0.07062697 0.1412539 0.9293730
[10,] 0.25042242 0.5008448 0.7495776
[11,] 0.20362283 0.4072457 0.7963772
[12,] 0.14649694 0.2929939 0.8535031
[13,] 0.12199956 0.2439991 0.8780004
[14,] 0.11473608 0.2294722 0.8852639
[15,] 0.10079042 0.2015808 0.8992096
[16,] 0.07155793 0.1431159 0.9284421
[17,] 0.09544892 0.1908978 0.9045511
[18,] 0.07364716 0.1472943 0.9263528
[19,] 0.06775412 0.1355082 0.9322459
[20,] 0.05961078 0.1192216 0.9403892
[21,] 0.04402225 0.0880445 0.9559777
[22,] 0.06998219 0.1399644 0.9300178
[23,] 0.05360958 0.1072192 0.9463904
[24,] 0.06253354 0.1250671 0.9374665
[25,] 0.05451405 0.1090281 0.9454859
[26,] 0.13822151 0.2764430 0.8617785
[27,] 0.12581538 0.2516308 0.8741846
[28,] 0.10546251 0.2109250 0.8945375
[29,] 0.09703076 0.1940615 0.9029692
[30,] 0.07595578 0.1519116 0.9240442
[31,] 0.12310303 0.2462061 0.8768970
[32,] 0.09847952 0.1969590 0.9015205
[33,] 0.08026817 0.1605363 0.9197318
[34,] 0.20182011 0.4036402 0.7981799
[35,] 0.17464105 0.3492821 0.8253589
[36,] 0.14473850 0.2894770 0.8552615
[37,] 0.17804428 0.3560886 0.8219557
[38,] 0.19112814 0.3822563 0.8088719
[39,] 0.21892319 0.4378464 0.7810768
[40,] 0.25872231 0.5174446 0.7412777
[41,] 0.23416928 0.4683386 0.7658307
[42,] 0.20539046 0.4107809 0.7946095
[43,] 0.28922552 0.5784510 0.7107745
[44,] 0.25648472 0.5129694 0.7435153
[45,] 0.23720174 0.4744035 0.7627983
[46,] 0.29237351 0.5847470 0.7076265
[47,] 0.26051590 0.5210318 0.7394841
[48,] 0.27621741 0.5524348 0.7237826
[49,] 0.48691474 0.9738295 0.5130853
[50,] 0.47481532 0.9496306 0.5251847
[51,] 0.49648992 0.9929798 0.5035101
[52,] 0.47887264 0.9577453 0.5211274
[53,] 0.45992285 0.9198457 0.5400771
[54,] 0.42617352 0.8523470 0.5738265
[55,] 0.42064309 0.8412862 0.5793569
[56,] 0.43016281 0.8603256 0.5698372
[57,] 0.39459279 0.7891856 0.6054072
[58,] 0.50098187 0.9980363 0.4990181
[59,] 0.55642264 0.8871547 0.4435774
[60,] 0.53302957 0.9339409 0.4669704
[61,] 0.51744507 0.9651099 0.4825549
[62,] 0.49698044 0.9939609 0.5030196
[63,] 0.50318697 0.9936261 0.4968130
[64,] 0.49699367 0.9939873 0.5030063
[65,] 0.47266776 0.9453355 0.5273322
[66,] 0.45103868 0.9020774 0.5489613
[67,] 0.47268036 0.9453607 0.5273196
[68,] 0.42523621 0.8504724 0.5747638
[69,] 0.41894038 0.8378808 0.5810596
[70,] 0.56858029 0.8628394 0.4314197
[71,] 0.61029570 0.7794086 0.3897043
[72,] 0.63001575 0.7399685 0.3699842
[73,] 0.59494494 0.8101101 0.4050551
[74,] 0.56242256 0.8751549 0.4375774
[75,] 0.51321214 0.9735757 0.4867879
[76,] 0.46455615 0.9291123 0.5354439
[77,] 0.41098961 0.8219792 0.5890104
[78,] 0.53464608 0.9307078 0.4653539
[79,] 0.46063663 0.9212733 0.5393634
[80,] 0.43167088 0.8633418 0.5683291
[81,] 0.45679790 0.9135958 0.5432021
[82,] 0.49799836 0.9959967 0.5020016
[83,] 0.42983064 0.8596613 0.5701694
[84,] 0.48276248 0.9655250 0.5172375
[85,] 0.37200158 0.7440032 0.6279984
[86,] 0.51439319 0.9712136 0.4856068
[87,] 0.36715665 0.7343133 0.6328434
> postscript(file="/var/www/html/freestat/rcomp/tmp/126es1291544884.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/www/html/freestat/rcomp/tmp/2dgdv1291544884.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/www/html/freestat/rcomp/tmp/3dgdv1291544884.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/www/html/freestat/rcomp/tmp/4dgdv1291544884.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/www/html/freestat/rcomp/tmp/567cy1291544884.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 = 96
Frequency = 1
1 2 3 4 5 6
200.602726 -316.845058 262.450357 -491.571015 228.288451 -539.823857
7 8 9 10 11 12
-130.514734 132.501870 -577.392251 157.720653 -283.239498 -272.933908
13 14 15 16 17 18
367.931561 -941.189049 145.260821 -99.155203 159.711117 -572.598816
19 20 21 22 23 24
192.585313 -79.150819 -751.207909 142.416937 -418.811638 -431.507709
25 26 27 28 29 30
140.820700 -734.824921 104.539207 -624.108288 -473.899637 -1063.199607
31 32 33 34 35 36
141.852076 -365.642198 -483.928076 13.614987 -786.121996 -225.536360
37 38 39 40 41 42
135.528819 -1004.663824 -274.427160 -135.874944 -666.954854 -536.349082
43 44 45 46 47 48
379.724613 -592.392251 -163.016926 -165.675234 -685.939102 98.111347
49 50 51 52 53 54
258.765483 -552.162058 134.846244 -375.792523 -786.829051 -67.756804
55 56 57 58 59 60
362.277041 162.797072 -2.411153 71.628483 -190.105180 476.762375
61 62 63 64 65 66
280.480243 -461.558158 669.520730 2.942374 8.620606 75.776722
67 68 69 70 71 72
-4.989722 481.446440 180.047019 481.368361 -181.207526 305.557684
73 74 75 76 77 78
598.421492 -393.982697 794.388285 -45.462028 545.825298 220.834025
79 80 81 82 83 84
381.751817 556.072350 560.119478 1041.962298 220.488758 59.794348
85 86 87 88 89 90
851.788941 -119.986018 696.054469 -33.973970 613.830279 -23.370453
91 92 93 94 95 96
754.275594 1115.117605 278.515366 1071.589869 266.404889 544.349671
> postscript(file="/var/www/html/freestat/rcomp/tmp/667cy1291544884.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 = 96
Frequency = 1
lag(myerror, k = 1) myerror
0 200.602726 NA
1 -316.845058 200.602726
2 262.450357 -316.845058
3 -491.571015 262.450357
4 228.288451 -491.571015
5 -539.823857 228.288451
6 -130.514734 -539.823857
7 132.501870 -130.514734
8 -577.392251 132.501870
9 157.720653 -577.392251
10 -283.239498 157.720653
11 -272.933908 -283.239498
12 367.931561 -272.933908
13 -941.189049 367.931561
14 145.260821 -941.189049
15 -99.155203 145.260821
16 159.711117 -99.155203
17 -572.598816 159.711117
18 192.585313 -572.598816
19 -79.150819 192.585313
20 -751.207909 -79.150819
21 142.416937 -751.207909
22 -418.811638 142.416937
23 -431.507709 -418.811638
24 140.820700 -431.507709
25 -734.824921 140.820700
26 104.539207 -734.824921
27 -624.108288 104.539207
28 -473.899637 -624.108288
29 -1063.199607 -473.899637
30 141.852076 -1063.199607
31 -365.642198 141.852076
32 -483.928076 -365.642198
33 13.614987 -483.928076
34 -786.121996 13.614987
35 -225.536360 -786.121996
36 135.528819 -225.536360
37 -1004.663824 135.528819
38 -274.427160 -1004.663824
39 -135.874944 -274.427160
40 -666.954854 -135.874944
41 -536.349082 -666.954854
42 379.724613 -536.349082
43 -592.392251 379.724613
44 -163.016926 -592.392251
45 -165.675234 -163.016926
46 -685.939102 -165.675234
47 98.111347 -685.939102
48 258.765483 98.111347
49 -552.162058 258.765483
50 134.846244 -552.162058
51 -375.792523 134.846244
52 -786.829051 -375.792523
53 -67.756804 -786.829051
54 362.277041 -67.756804
55 162.797072 362.277041
56 -2.411153 162.797072
57 71.628483 -2.411153
58 -190.105180 71.628483
59 476.762375 -190.105180
60 280.480243 476.762375
61 -461.558158 280.480243
62 669.520730 -461.558158
63 2.942374 669.520730
64 8.620606 2.942374
65 75.776722 8.620606
66 -4.989722 75.776722
67 481.446440 -4.989722
68 180.047019 481.446440
69 481.368361 180.047019
70 -181.207526 481.368361
71 305.557684 -181.207526
72 598.421492 305.557684
73 -393.982697 598.421492
74 794.388285 -393.982697
75 -45.462028 794.388285
76 545.825298 -45.462028
77 220.834025 545.825298
78 381.751817 220.834025
79 556.072350 381.751817
80 560.119478 556.072350
81 1041.962298 560.119478
82 220.488758 1041.962298
83 59.794348 220.488758
84 851.788941 59.794348
85 -119.986018 851.788941
86 696.054469 -119.986018
87 -33.973970 696.054469
88 613.830279 -33.973970
89 -23.370453 613.830279
90 754.275594 -23.370453
91 1115.117605 754.275594
92 278.515366 1115.117605
93 1071.589869 278.515366
94 266.404889 1071.589869
95 544.349671 266.404889
96 NA 544.349671
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -316.845058 200.602726
[2,] 262.450357 -316.845058
[3,] -491.571015 262.450357
[4,] 228.288451 -491.571015
[5,] -539.823857 228.288451
[6,] -130.514734 -539.823857
[7,] 132.501870 -130.514734
[8,] -577.392251 132.501870
[9,] 157.720653 -577.392251
[10,] -283.239498 157.720653
[11,] -272.933908 -283.239498
[12,] 367.931561 -272.933908
[13,] -941.189049 367.931561
[14,] 145.260821 -941.189049
[15,] -99.155203 145.260821
[16,] 159.711117 -99.155203
[17,] -572.598816 159.711117
[18,] 192.585313 -572.598816
[19,] -79.150819 192.585313
[20,] -751.207909 -79.150819
[21,] 142.416937 -751.207909
[22,] -418.811638 142.416937
[23,] -431.507709 -418.811638
[24,] 140.820700 -431.507709
[25,] -734.824921 140.820700
[26,] 104.539207 -734.824921
[27,] -624.108288 104.539207
[28,] -473.899637 -624.108288
[29,] -1063.199607 -473.899637
[30,] 141.852076 -1063.199607
[31,] -365.642198 141.852076
[32,] -483.928076 -365.642198
[33,] 13.614987 -483.928076
[34,] -786.121996 13.614987
[35,] -225.536360 -786.121996
[36,] 135.528819 -225.536360
[37,] -1004.663824 135.528819
[38,] -274.427160 -1004.663824
[39,] -135.874944 -274.427160
[40,] -666.954854 -135.874944
[41,] -536.349082 -666.954854
[42,] 379.724613 -536.349082
[43,] -592.392251 379.724613
[44,] -163.016926 -592.392251
[45,] -165.675234 -163.016926
[46,] -685.939102 -165.675234
[47,] 98.111347 -685.939102
[48,] 258.765483 98.111347
[49,] -552.162058 258.765483
[50,] 134.846244 -552.162058
[51,] -375.792523 134.846244
[52,] -786.829051 -375.792523
[53,] -67.756804 -786.829051
[54,] 362.277041 -67.756804
[55,] 162.797072 362.277041
[56,] -2.411153 162.797072
[57,] 71.628483 -2.411153
[58,] -190.105180 71.628483
[59,] 476.762375 -190.105180
[60,] 280.480243 476.762375
[61,] -461.558158 280.480243
[62,] 669.520730 -461.558158
[63,] 2.942374 669.520730
[64,] 8.620606 2.942374
[65,] 75.776722 8.620606
[66,] -4.989722 75.776722
[67,] 481.446440 -4.989722
[68,] 180.047019 481.446440
[69,] 481.368361 180.047019
[70,] -181.207526 481.368361
[71,] 305.557684 -181.207526
[72,] 598.421492 305.557684
[73,] -393.982697 598.421492
[74,] 794.388285 -393.982697
[75,] -45.462028 794.388285
[76,] 545.825298 -45.462028
[77,] 220.834025 545.825298
[78,] 381.751817 220.834025
[79,] 556.072350 381.751817
[80,] 560.119478 556.072350
[81,] 1041.962298 560.119478
[82,] 220.488758 1041.962298
[83,] 59.794348 220.488758
[84,] 851.788941 59.794348
[85,] -119.986018 851.788941
[86,] 696.054469 -119.986018
[87,] -33.973970 696.054469
[88,] 613.830279 -33.973970
[89,] -23.370453 613.830279
[90,] 754.275594 -23.370453
[91,] 1115.117605 754.275594
[92,] 278.515366 1115.117605
[93,] 1071.589869 278.515366
[94,] 266.404889 1071.589869
[95,] 544.349671 266.404889
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -316.845058 200.602726
2 262.450357 -316.845058
3 -491.571015 262.450357
4 228.288451 -491.571015
5 -539.823857 228.288451
6 -130.514734 -539.823857
7 132.501870 -130.514734
8 -577.392251 132.501870
9 157.720653 -577.392251
10 -283.239498 157.720653
11 -272.933908 -283.239498
12 367.931561 -272.933908
13 -941.189049 367.931561
14 145.260821 -941.189049
15 -99.155203 145.260821
16 159.711117 -99.155203
17 -572.598816 159.711117
18 192.585313 -572.598816
19 -79.150819 192.585313
20 -751.207909 -79.150819
21 142.416937 -751.207909
22 -418.811638 142.416937
23 -431.507709 -418.811638
24 140.820700 -431.507709
25 -734.824921 140.820700
26 104.539207 -734.824921
27 -624.108288 104.539207
28 -473.899637 -624.108288
29 -1063.199607 -473.899637
30 141.852076 -1063.199607
31 -365.642198 141.852076
32 -483.928076 -365.642198
33 13.614987 -483.928076
34 -786.121996 13.614987
35 -225.536360 -786.121996
36 135.528819 -225.536360
37 -1004.663824 135.528819
38 -274.427160 -1004.663824
39 -135.874944 -274.427160
40 -666.954854 -135.874944
41 -536.349082 -666.954854
42 379.724613 -536.349082
43 -592.392251 379.724613
44 -163.016926 -592.392251
45 -165.675234 -163.016926
46 -685.939102 -165.675234
47 98.111347 -685.939102
48 258.765483 98.111347
49 -552.162058 258.765483
50 134.846244 -552.162058
51 -375.792523 134.846244
52 -786.829051 -375.792523
53 -67.756804 -786.829051
54 362.277041 -67.756804
55 162.797072 362.277041
56 -2.411153 162.797072
57 71.628483 -2.411153
58 -190.105180 71.628483
59 476.762375 -190.105180
60 280.480243 476.762375
61 -461.558158 280.480243
62 669.520730 -461.558158
63 2.942374 669.520730
64 8.620606 2.942374
65 75.776722 8.620606
66 -4.989722 75.776722
67 481.446440 -4.989722
68 180.047019 481.446440
69 481.368361 180.047019
70 -181.207526 481.368361
71 305.557684 -181.207526
72 598.421492 305.557684
73 -393.982697 598.421492
74 794.388285 -393.982697
75 -45.462028 794.388285
76 545.825298 -45.462028
77 220.834025 545.825298
78 381.751817 220.834025
79 556.072350 381.751817
80 560.119478 556.072350
81 1041.962298 560.119478
82 220.488758 1041.962298
83 59.794348 220.488758
84 851.788941 59.794348
85 -119.986018 851.788941
86 696.054469 -119.986018
87 -33.973970 696.054469
88 613.830279 -33.973970
89 -23.370453 613.830279
90 754.275594 -23.370453
91 1115.117605 754.275594
92 278.515366 1115.117605
93 1071.589869 278.515366
94 266.404889 1071.589869
95 544.349671 266.404889
> 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/freestat/rcomp/tmp/7gguj1291544884.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/www/html/freestat/rcomp/tmp/8gguj1291544884.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/www/html/freestat/rcomp/tmp/998bm1291544884.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/www/html/freestat/rcomp/tmp/1098bm1291544884.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/www/html/freestat/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/html/freestat/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/freestat/rcomp/tmp/11u8rs1291544884.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/freestat/rcomp/tmp/12y88g1291544884.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/freestat/rcomp/tmp/13c05p1291544884.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/freestat/rcomp/tmp/14fjmd1291544884.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/freestat/rcomp/tmp/15jjlj1291544884.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/freestat/rcomp/tmp/16ftir1291544884.tab")
+ }
>
> try(system("convert tmp/126es1291544884.ps tmp/126es1291544884.png",intern=TRUE))
character(0)
> try(system("convert tmp/2dgdv1291544884.ps tmp/2dgdv1291544884.png",intern=TRUE))
character(0)
> try(system("convert tmp/3dgdv1291544884.ps tmp/3dgdv1291544884.png",intern=TRUE))
character(0)
> try(system("convert tmp/4dgdv1291544884.ps tmp/4dgdv1291544884.png",intern=TRUE))
character(0)
> try(system("convert tmp/567cy1291544884.ps tmp/567cy1291544884.png",intern=TRUE))
character(0)
> try(system("convert tmp/667cy1291544884.ps tmp/667cy1291544884.png",intern=TRUE))
character(0)
> try(system("convert tmp/7gguj1291544884.ps tmp/7gguj1291544884.png",intern=TRUE))
character(0)
> try(system("convert tmp/8gguj1291544884.ps tmp/8gguj1291544884.png",intern=TRUE))
character(0)
> try(system("convert tmp/998bm1291544884.ps tmp/998bm1291544884.png",intern=TRUE))
character(0)
> try(system("convert tmp/1098bm1291544884.ps tmp/1098bm1291544884.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
4.326 2.559 4.619