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(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 = 'Linear Trend'
> par2 = 'Include Monthly 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 M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 9769 1579 1 0 0 0 0 0 0 0 0 0 0 1
2 9321 2146 0 1 0 0 0 0 0 0 0 0 0 2
3 9939 2462 0 0 1 0 0 0 0 0 0 0 0 3
4 9336 3695 0 0 0 1 0 0 0 0 0 0 0 4
5 10195 4831 0 0 0 0 1 0 0 0 0 0 0 5
6 9464 5134 0 0 0 0 0 1 0 0 0 0 0 6
7 10010 6250 0 0 0 0 0 0 1 0 0 0 0 7
8 10213 5760 0 0 0 0 0 0 0 1 0 0 0 8
9 9563 6249 0 0 0 0 0 0 0 0 1 0 0 9
10 9890 2917 0 0 0 0 0 0 0 0 0 1 0 10
11 9305 1741 0 0 0 0 0 0 0 0 0 0 1 11
12 9391 2359 0 0 0 0 0 0 0 0 0 0 0 12
13 9928 1511 1 0 0 0 0 0 0 0 0 0 0 13
14 8686 2059 0 1 0 0 0 0 0 0 0 0 0 14
15 9843 2635 0 0 1 0 0 0 0 0 0 0 0 15
16 9627 2867 0 0 0 1 0 0 0 0 0 0 0 16
17 10074 4403 0 0 0 0 1 0 0 0 0 0 0 17
18 9503 5720 0 0 0 0 0 1 0 0 0 0 0 18
19 10119 4502 0 0 0 0 0 0 1 0 0 0 0 19
20 10000 5749 0 0 0 0 0 0 0 1 0 0 0 20
21 9313 5627 0 0 0 0 0 0 0 0 1 0 0 21
22 9866 2846 0 0 0 0 0 0 0 0 0 1 0 22
23 9172 1762 0 0 0 0 0 0 0 0 0 0 1 23
24 9241 2429 0 0 0 0 0 0 0 0 0 0 0 24
25 9659 1169 1 0 0 0 0 0 0 0 0 0 0 25
26 8904 2154 0 1 0 0 0 0 0 0 0 0 0 26
27 9755 2249 0 0 1 0 0 0 0 0 0 0 0 27
28 9080 2687 0 0 0 1 0 0 0 0 0 0 0 28
29 9435 4359 0 0 0 0 1 0 0 0 0 0 0 29
30 8971 5382 0 0 0 0 0 1 0 0 0 0 0 30
31 10063 4459 0 0 0 0 0 0 1 0 0 0 0 31
32 9793 6398 0 0 0 0 0 0 0 1 0 0 0 32
33 9454 4596 0 0 0 0 0 0 0 0 1 0 0 33
34 9759 3024 0 0 0 0 0 0 0 0 0 1 0 34
35 8820 1887 0 0 0 0 0 0 0 0 0 0 1 35
36 9403 2070 0 0 0 0 0 0 0 0 0 0 0 36
37 9676 1351 1 0 0 0 0 0 0 0 0 0 0 37
38 8642 2218 0 1 0 0 0 0 0 0 0 0 0 38
39 9402 2461 0 0 1 0 0 0 0 0 0 0 0 39
40 9610 3028 0 0 0 1 0 0 0 0 0 0 0 40
41 9294 4784 0 0 0 0 1 0 0 0 0 0 0 41
42 9448 4975 0 0 0 0 0 1 0 0 0 0 0 42
43 10319 4607 0 0 0 0 0 0 1 0 0 0 0 43
44 9548 6249 0 0 0 0 0 0 0 1 0 0 0 44
45 9801 4809 0 0 0 0 0 0 0 0 1 0 0 45
46 9596 3157 0 0 0 0 0 0 0 0 0 1 0 46
47 8923 1910 0 0 0 0 0 0 0 0 0 0 1 47
48 9746 2228 0 0 0 0 0 0 0 0 0 0 0 48
49 9829 1594 1 0 0 0 0 0 0 0 0 0 0 49
50 9125 2467 0 1 0 0 0 0 0 0 0 0 0 50
51 9782 2222 0 0 1 0 0 0 0 0 0 0 0 51
52 9441 3607 0 0 0 1 0 0 0 0 0 0 0 52
53 9162 4685 0 0 0 0 1 0 0 0 0 0 0 53
54 9915 4962 0 0 0 0 0 1 0 0 0 0 0 54
55 10444 5770 0 0 0 0 0 0 1 0 0 0 0 55
56 10209 5480 0 0 0 0 0 0 0 1 0 0 0 56
57 9985 5000 0 0 0 0 0 0 0 0 1 0 0 57
58 9842 3228 0 0 0 0 0 0 0 0 0 1 0 58
59 9429 1993 0 0 0 0 0 0 0 0 0 0 1 59
60 10132 2288 0 0 0 0 0 0 0 0 0 0 0 60
61 9849 1580 1 0 0 0 0 0 0 0 0 0 0 61
62 9172 2111 0 1 0 0 0 0 0 0 0 0 0 62
63 10313 2192 0 0 1 0 0 0 0 0 0 0 0 63
64 9819 3601 0 0 0 1 0 0 0 0 0 0 0 64
65 9955 4665 0 0 0 0 1 0 0 0 0 0 0 65
66 10048 4876 0 0 0 0 0 1 0 0 0 0 0 66
67 10082 5813 0 0 0 0 0 0 1 0 0 0 0 67
68 10541 5589 0 0 0 0 0 0 0 1 0 0 0 68
69 10208 5331 0 0 0 0 0 0 0 0 1 0 0 69
70 10233 3075 0 0 0 0 0 0 0 0 0 1 0 70
71 9439 2002 0 0 0 0 0 0 0 0 0 0 1 71
72 9963 2306 0 0 0 0 0 0 0 0 0 0 0 72
73 10158 1507 1 0 0 0 0 0 0 0 0 0 0 73
74 9225 1992 0 1 0 0 0 0 0 0 0 0 0 74
75 10474 2487 0 0 1 0 0 0 0 0 0 0 0 75
76 9757 3490 0 0 0 1 0 0 0 0 0 0 0 76
77 10490 4647 0 0 0 0 1 0 0 0 0 0 0 77
78 10281 5594 0 0 0 0 0 1 0 0 0 0 0 78
79 10444 5611 0 0 0 0 0 0 1 0 0 0 0 79
80 10640 5788 0 0 0 0 0 0 0 1 0 0 0 80
81 10695 6204 0 0 0 0 0 0 0 0 1 0 0 81
82 10786 3013 0 0 0 0 0 0 0 0 0 1 0 82
83 9832 1931 0 0 0 0 0 0 0 0 0 0 1 83
84 9747 2549 0 0 0 0 0 0 0 0 0 0 0 84
85 10411 1504 1 0 0 0 0 0 0 0 0 0 0 85
86 9511 2090 0 1 0 0 0 0 0 0 0 0 0 86
87 10402 2702 0 0 1 0 0 0 0 0 0 0 0 87
88 9701 2939 0 0 0 1 0 0 0 0 0 0 0 88
89 10540 4500 0 0 0 0 1 0 0 0 0 0 0 89
90 10112 6208 0 0 0 0 0 1 0 0 0 0 0 90
91 10915 6415 0 0 0 0 0 0 1 0 0 0 0 91
92 11183 5657 0 0 0 0 0 0 0 1 0 0 0 92
93 10384 5964 0 0 0 0 0 0 0 0 1 0 0 93
94 10834 3163 0 0 0 0 0 0 0 0 0 1 0 94
95 9886 1997 0 0 0 0 0 0 0 0 0 0 1 95
96 10216 2422 0 0 0 0 0 0 0 0 0 0 0 96
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) huwelijken M1 M2 M3 M4
9198.18673 0.01321 293.35422 -561.53253 341.10358 -121.28678
M5 M6 M7 M8 M9 M10
198.08998 3.56981 575.09294 526.83724 181.83353 379.89563
M11 t
-364.18882 9.27576
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-787.78 -168.02 20.72 150.20 688.53
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 9198.18673 229.71239 40.042 < 2e-16 ***
huwelijken 0.01321 0.08823 0.150 0.881347
M1 293.35422 166.28689 1.764 0.081432 .
M2 -561.53253 149.69917 -3.751 0.000327 ***
M3 341.10358 149.33640 2.284 0.024949 *
M4 -121.28678 169.82913 -0.714 0.477150
M5 198.08998 251.17373 0.789 0.432587
M6 3.56981 306.56970 0.012 0.990738
M7 575.09294 311.94906 1.844 0.068862 .
M8 526.83724 343.61966 1.533 0.129076
M9 181.83353 315.00813 0.577 0.565363
M10 379.89563 161.88425 2.347 0.021353 *
M11 -364.18882 153.35203 -2.375 0.019891 *
t 9.27576 1.12011 8.281 1.94e-12 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 297.3 on 82 degrees of freedom
Multiple R-squared: 0.7099, Adjusted R-squared: 0.6639
F-statistic: 15.44 on 13 and 82 DF, p-value: < 2.2e-16
> 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.67417404 0.65165193 0.32582596
[2,] 0.65159093 0.69681814 0.34840907
[3,] 0.52096052 0.95807896 0.47903948
[4,] 0.42409746 0.84819492 0.57590254
[5,] 0.35681008 0.71362016 0.64318992
[6,] 0.27661797 0.55323594 0.72338203
[7,] 0.22071744 0.44143488 0.77928256
[8,] 0.15123616 0.30247232 0.84876384
[9,] 0.10595630 0.21191260 0.89404370
[10,] 0.08333737 0.16667474 0.91666263
[11,] 0.05808497 0.11616993 0.94191503
[12,] 0.06763804 0.13527608 0.93236196
[13,] 0.19232214 0.38464428 0.80767786
[14,] 0.19853863 0.39707726 0.80146137
[15,] 0.16548188 0.33096376 0.83451812
[16,] 0.12348495 0.24696990 0.87651505
[17,] 0.09927126 0.19854252 0.90072874
[18,] 0.08172393 0.16344786 0.91827607
[19,] 0.06178610 0.12357221 0.93821390
[20,] 0.06679306 0.13358611 0.93320694
[21,] 0.05878193 0.11756386 0.94121807
[22,] 0.03970429 0.07940859 0.96029571
[23,] 0.03300849 0.06601698 0.96699151
[24,] 0.13017981 0.26035962 0.86982019
[25,] 0.13427520 0.26855039 0.86572480
[26,] 0.15479392 0.30958784 0.84520608
[27,] 0.21163359 0.42326717 0.78836641
[28,] 0.25298616 0.50597232 0.74701384
[29,] 0.31330260 0.62660520 0.68669740
[30,] 0.29613237 0.59226474 0.70386763
[31,] 0.26514497 0.53028993 0.73485503
[32,] 0.40552593 0.81105187 0.59447407
[33,] 0.40022448 0.80044897 0.59977552
[34,] 0.47225057 0.94450113 0.52774943
[35,] 0.44163648 0.88327295 0.55836352
[36,] 0.40647628 0.81295256 0.59352372
[37,] 0.74526471 0.50947058 0.25473529
[38,] 0.82395758 0.35208484 0.17604242
[39,] 0.87767641 0.24464717 0.12232359
[40,] 0.87568741 0.24862518 0.12431259
[41,] 0.87563937 0.24872126 0.12436063
[42,] 0.91148268 0.17703465 0.08851732
[43,] 0.90053754 0.19892491 0.09946246
[44,] 0.97088094 0.05823813 0.02911906
[45,] 0.96061129 0.07877742 0.03938871
[46,] 0.94431012 0.11137975 0.05568988
[47,] 0.94308411 0.11383178 0.05691589
[48,] 0.96075022 0.07849956 0.03924978
[49,] 0.95648069 0.08703861 0.04351931
[50,] 0.95153680 0.09692640 0.04846320
[51,] 0.95148644 0.09702712 0.04851356
[52,] 0.93717214 0.12565572 0.06282786
[53,] 0.91589809 0.16820381 0.08410191
[54,] 0.92987228 0.14025544 0.07012772
[55,] 0.91907692 0.16184616 0.08092308
[56,] 0.89027777 0.21944446 0.10972223
[57,] 0.84178963 0.31642073 0.15821037
[58,] 0.78138973 0.43722055 0.21861027
[59,] 0.73168874 0.53662251 0.26831126
[60,] 0.62883822 0.74232356 0.37116178
[61,] 0.53228810 0.93542381 0.46771190
[62,] 0.57250625 0.85498750 0.42749375
[63,] 0.41215548 0.82431095 0.58784452
> postscript(file="/var/www/html/rcomp/tmp/15dmw1292080178.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/rcomp/tmp/2gm4z1292080178.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/rcomp/tmp/3gm4z1292080178.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/rcomp/tmp/4gm4z1292080178.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/rcomp/tmp/5rdl21292080178.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
247.324951 637.445965 339.359779 173.186646 688.527749 138.769569
7 8 9 10 11 12
89.228496 337.681256 16.949590 180.626917 345.970373 50.342109
13 14 15 16 17 18
295.914069 -107.713930 129.765325 363.815242 461.872408 58.719451
19 20 21 22 23 24
110.010144 13.517413 -336.143041 46.255664 101.383815 -211.891732
25 26 27 28 29 30
-79.877318 -2.278016 -64.444829 -292.116139 -287.855510 -580.124775
31 32 33 34 35 36
-56.730985 -313.364924 -292.832849 -174.404839 -363.576566 -156.458551
37 38 39 40 41 42
-176.590660 -376.432598 -531.554466 122.070154 -545.778844 -209.057522
43 44 45 46 47 48
86.004807 -667.705810 -59.955696 -450.470899 -372.189544 73.145143
49 50 51 52 53 54
-138.109803 -8.030999 -259.706466 -165.887494 -787.780221 146.805055
55 56 57 58 59 60
84.332613 -107.856595 10.212073 -316.717949 21.404887 347.043401
61 62 63 64 65 66
-229.234017 -67.637448 160.380677 100.882613 -105.825176 169.631949
67 68 69 70 71 72
-389.544562 111.394381 117.530465 -35.005995 -80.023153 66.496472
73 74 75 76 77 78
-30.578850 -124.374629 206.174623 -70.960246 318.103450 281.838133
79 80 81 82 83 84
-136.185326 96.456471 481.689124 407.503863 202.605594 -264.022670
85 86 87 88 89 90
111.151628 49.021655 20.025356 -230.990776 258.736144 -6.581860
91 92 93 94 95 96
212.884812 529.877809 62.550333 342.213236 144.424593 95.345828
> postscript(file="/var/www/html/rcomp/tmp/6rdl21292080178.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 247.324951 NA
1 637.445965 247.324951
2 339.359779 637.445965
3 173.186646 339.359779
4 688.527749 173.186646
5 138.769569 688.527749
6 89.228496 138.769569
7 337.681256 89.228496
8 16.949590 337.681256
9 180.626917 16.949590
10 345.970373 180.626917
11 50.342109 345.970373
12 295.914069 50.342109
13 -107.713930 295.914069
14 129.765325 -107.713930
15 363.815242 129.765325
16 461.872408 363.815242
17 58.719451 461.872408
18 110.010144 58.719451
19 13.517413 110.010144
20 -336.143041 13.517413
21 46.255664 -336.143041
22 101.383815 46.255664
23 -211.891732 101.383815
24 -79.877318 -211.891732
25 -2.278016 -79.877318
26 -64.444829 -2.278016
27 -292.116139 -64.444829
28 -287.855510 -292.116139
29 -580.124775 -287.855510
30 -56.730985 -580.124775
31 -313.364924 -56.730985
32 -292.832849 -313.364924
33 -174.404839 -292.832849
34 -363.576566 -174.404839
35 -156.458551 -363.576566
36 -176.590660 -156.458551
37 -376.432598 -176.590660
38 -531.554466 -376.432598
39 122.070154 -531.554466
40 -545.778844 122.070154
41 -209.057522 -545.778844
42 86.004807 -209.057522
43 -667.705810 86.004807
44 -59.955696 -667.705810
45 -450.470899 -59.955696
46 -372.189544 -450.470899
47 73.145143 -372.189544
48 -138.109803 73.145143
49 -8.030999 -138.109803
50 -259.706466 -8.030999
51 -165.887494 -259.706466
52 -787.780221 -165.887494
53 146.805055 -787.780221
54 84.332613 146.805055
55 -107.856595 84.332613
56 10.212073 -107.856595
57 -316.717949 10.212073
58 21.404887 -316.717949
59 347.043401 21.404887
60 -229.234017 347.043401
61 -67.637448 -229.234017
62 160.380677 -67.637448
63 100.882613 160.380677
64 -105.825176 100.882613
65 169.631949 -105.825176
66 -389.544562 169.631949
67 111.394381 -389.544562
68 117.530465 111.394381
69 -35.005995 117.530465
70 -80.023153 -35.005995
71 66.496472 -80.023153
72 -30.578850 66.496472
73 -124.374629 -30.578850
74 206.174623 -124.374629
75 -70.960246 206.174623
76 318.103450 -70.960246
77 281.838133 318.103450
78 -136.185326 281.838133
79 96.456471 -136.185326
80 481.689124 96.456471
81 407.503863 481.689124
82 202.605594 407.503863
83 -264.022670 202.605594
84 111.151628 -264.022670
85 49.021655 111.151628
86 20.025356 49.021655
87 -230.990776 20.025356
88 258.736144 -230.990776
89 -6.581860 258.736144
90 212.884812 -6.581860
91 529.877809 212.884812
92 62.550333 529.877809
93 342.213236 62.550333
94 144.424593 342.213236
95 95.345828 144.424593
96 NA 95.345828
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 637.445965 247.324951
[2,] 339.359779 637.445965
[3,] 173.186646 339.359779
[4,] 688.527749 173.186646
[5,] 138.769569 688.527749
[6,] 89.228496 138.769569
[7,] 337.681256 89.228496
[8,] 16.949590 337.681256
[9,] 180.626917 16.949590
[10,] 345.970373 180.626917
[11,] 50.342109 345.970373
[12,] 295.914069 50.342109
[13,] -107.713930 295.914069
[14,] 129.765325 -107.713930
[15,] 363.815242 129.765325
[16,] 461.872408 363.815242
[17,] 58.719451 461.872408
[18,] 110.010144 58.719451
[19,] 13.517413 110.010144
[20,] -336.143041 13.517413
[21,] 46.255664 -336.143041
[22,] 101.383815 46.255664
[23,] -211.891732 101.383815
[24,] -79.877318 -211.891732
[25,] -2.278016 -79.877318
[26,] -64.444829 -2.278016
[27,] -292.116139 -64.444829
[28,] -287.855510 -292.116139
[29,] -580.124775 -287.855510
[30,] -56.730985 -580.124775
[31,] -313.364924 -56.730985
[32,] -292.832849 -313.364924
[33,] -174.404839 -292.832849
[34,] -363.576566 -174.404839
[35,] -156.458551 -363.576566
[36,] -176.590660 -156.458551
[37,] -376.432598 -176.590660
[38,] -531.554466 -376.432598
[39,] 122.070154 -531.554466
[40,] -545.778844 122.070154
[41,] -209.057522 -545.778844
[42,] 86.004807 -209.057522
[43,] -667.705810 86.004807
[44,] -59.955696 -667.705810
[45,] -450.470899 -59.955696
[46,] -372.189544 -450.470899
[47,] 73.145143 -372.189544
[48,] -138.109803 73.145143
[49,] -8.030999 -138.109803
[50,] -259.706466 -8.030999
[51,] -165.887494 -259.706466
[52,] -787.780221 -165.887494
[53,] 146.805055 -787.780221
[54,] 84.332613 146.805055
[55,] -107.856595 84.332613
[56,] 10.212073 -107.856595
[57,] -316.717949 10.212073
[58,] 21.404887 -316.717949
[59,] 347.043401 21.404887
[60,] -229.234017 347.043401
[61,] -67.637448 -229.234017
[62,] 160.380677 -67.637448
[63,] 100.882613 160.380677
[64,] -105.825176 100.882613
[65,] 169.631949 -105.825176
[66,] -389.544562 169.631949
[67,] 111.394381 -389.544562
[68,] 117.530465 111.394381
[69,] -35.005995 117.530465
[70,] -80.023153 -35.005995
[71,] 66.496472 -80.023153
[72,] -30.578850 66.496472
[73,] -124.374629 -30.578850
[74,] 206.174623 -124.374629
[75,] -70.960246 206.174623
[76,] 318.103450 -70.960246
[77,] 281.838133 318.103450
[78,] -136.185326 281.838133
[79,] 96.456471 -136.185326
[80,] 481.689124 96.456471
[81,] 407.503863 481.689124
[82,] 202.605594 407.503863
[83,] -264.022670 202.605594
[84,] 111.151628 -264.022670
[85,] 49.021655 111.151628
[86,] 20.025356 49.021655
[87,] -230.990776 20.025356
[88,] 258.736144 -230.990776
[89,] -6.581860 258.736144
[90,] 212.884812 -6.581860
[91,] 529.877809 212.884812
[92,] 62.550333 529.877809
[93,] 342.213236 62.550333
[94,] 144.424593 342.213236
[95,] 95.345828 144.424593
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 637.445965 247.324951
2 339.359779 637.445965
3 173.186646 339.359779
4 688.527749 173.186646
5 138.769569 688.527749
6 89.228496 138.769569
7 337.681256 89.228496
8 16.949590 337.681256
9 180.626917 16.949590
10 345.970373 180.626917
11 50.342109 345.970373
12 295.914069 50.342109
13 -107.713930 295.914069
14 129.765325 -107.713930
15 363.815242 129.765325
16 461.872408 363.815242
17 58.719451 461.872408
18 110.010144 58.719451
19 13.517413 110.010144
20 -336.143041 13.517413
21 46.255664 -336.143041
22 101.383815 46.255664
23 -211.891732 101.383815
24 -79.877318 -211.891732
25 -2.278016 -79.877318
26 -64.444829 -2.278016
27 -292.116139 -64.444829
28 -287.855510 -292.116139
29 -580.124775 -287.855510
30 -56.730985 -580.124775
31 -313.364924 -56.730985
32 -292.832849 -313.364924
33 -174.404839 -292.832849
34 -363.576566 -174.404839
35 -156.458551 -363.576566
36 -176.590660 -156.458551
37 -376.432598 -176.590660
38 -531.554466 -376.432598
39 122.070154 -531.554466
40 -545.778844 122.070154
41 -209.057522 -545.778844
42 86.004807 -209.057522
43 -667.705810 86.004807
44 -59.955696 -667.705810
45 -450.470899 -59.955696
46 -372.189544 -450.470899
47 73.145143 -372.189544
48 -138.109803 73.145143
49 -8.030999 -138.109803
50 -259.706466 -8.030999
51 -165.887494 -259.706466
52 -787.780221 -165.887494
53 146.805055 -787.780221
54 84.332613 146.805055
55 -107.856595 84.332613
56 10.212073 -107.856595
57 -316.717949 10.212073
58 21.404887 -316.717949
59 347.043401 21.404887
60 -229.234017 347.043401
61 -67.637448 -229.234017
62 160.380677 -67.637448
63 100.882613 160.380677
64 -105.825176 100.882613
65 169.631949 -105.825176
66 -389.544562 169.631949
67 111.394381 -389.544562
68 117.530465 111.394381
69 -35.005995 117.530465
70 -80.023153 -35.005995
71 66.496472 -80.023153
72 -30.578850 66.496472
73 -124.374629 -30.578850
74 206.174623 -124.374629
75 -70.960246 206.174623
76 318.103450 -70.960246
77 281.838133 318.103450
78 -136.185326 281.838133
79 96.456471 -136.185326
80 481.689124 96.456471
81 407.503863 481.689124
82 202.605594 407.503863
83 -264.022670 202.605594
84 111.151628 -264.022670
85 49.021655 111.151628
86 20.025356 49.021655
87 -230.990776 20.025356
88 258.736144 -230.990776
89 -6.581860 258.736144
90 212.884812 -6.581860
91 529.877809 212.884812
92 62.550333 529.877809
93 342.213236 62.550333
94 144.424593 342.213236
95 95.345828 144.424593
> 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/724k51292080178.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/rcomp/tmp/824k51292080178.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/rcomp/tmp/9uwk81292080178.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/rcomp/tmp/10uwk81292080178.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/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/11gw0w1292080178.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/12jfy21292080178.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/13qgvv1292080178.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/14tguj1292080178.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/15mptm1292080178.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/16ihrd1292080178.tab")
+ }
>
> try(system("convert tmp/15dmw1292080178.ps tmp/15dmw1292080178.png",intern=TRUE))
character(0)
> try(system("convert tmp/2gm4z1292080178.ps tmp/2gm4z1292080178.png",intern=TRUE))
character(0)
> try(system("convert tmp/3gm4z1292080178.ps tmp/3gm4z1292080178.png",intern=TRUE))
character(0)
> try(system("convert tmp/4gm4z1292080178.ps tmp/4gm4z1292080178.png",intern=TRUE))
character(0)
> try(system("convert tmp/5rdl21292080178.ps tmp/5rdl21292080178.png",intern=TRUE))
character(0)
> try(system("convert tmp/6rdl21292080178.ps tmp/6rdl21292080178.png",intern=TRUE))
character(0)
> try(system("convert tmp/724k51292080178.ps tmp/724k51292080178.png",intern=TRUE))
character(0)
> try(system("convert tmp/824k51292080178.ps tmp/824k51292080178.png",intern=TRUE))
character(0)
> try(system("convert tmp/9uwk81292080178.ps tmp/9uwk81292080178.png",intern=TRUE))
character(0)
> try(system("convert tmp/10uwk81292080178.ps tmp/10uwk81292080178.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
2.943 1.644 6.958