R version 2.12.0 (2010-10-15)
Copyright (C) 2010 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: i486-pc-linux-gnu (32-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> x <- array(list(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 = '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
> 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
1 9769 1579 1 0 0 0 0 0 0 0 0 0 0
2 9321 2146 0 1 0 0 0 0 0 0 0 0 0
3 9939 2462 0 0 1 0 0 0 0 0 0 0 0
4 9336 3695 0 0 0 1 0 0 0 0 0 0 0
5 10195 4831 0 0 0 0 1 0 0 0 0 0 0
6 9464 5134 0 0 0 0 0 1 0 0 0 0 0
7 10010 6250 0 0 0 0 0 0 1 0 0 0 0
8 10213 5760 0 0 0 0 0 0 0 1 0 0 0
9 9563 6249 0 0 0 0 0 0 0 0 1 0 0
10 9890 2917 0 0 0 0 0 0 0 0 0 1 0
11 9305 1741 0 0 0 0 0 0 0 0 0 0 1
12 9391 2359 0 0 0 0 0 0 0 0 0 0 0
13 9928 1511 1 0 0 0 0 0 0 0 0 0 0
14 8686 2059 0 1 0 0 0 0 0 0 0 0 0
15 9843 2635 0 0 1 0 0 0 0 0 0 0 0
16 9627 2867 0 0 0 1 0 0 0 0 0 0 0
17 10074 4403 0 0 0 0 1 0 0 0 0 0 0
18 9503 5720 0 0 0 0 0 1 0 0 0 0 0
19 10119 4502 0 0 0 0 0 0 1 0 0 0 0
20 10000 5749 0 0 0 0 0 0 0 1 0 0 0
21 9313 5627 0 0 0 0 0 0 0 0 1 0 0
22 9866 2846 0 0 0 0 0 0 0 0 0 1 0
23 9172 1762 0 0 0 0 0 0 0 0 0 0 1
24 9241 2429 0 0 0 0 0 0 0 0 0 0 0
25 9659 1169 1 0 0 0 0 0 0 0 0 0 0
26 8904 2154 0 1 0 0 0 0 0 0 0 0 0
27 9755 2249 0 0 1 0 0 0 0 0 0 0 0
28 9080 2687 0 0 0 1 0 0 0 0 0 0 0
29 9435 4359 0 0 0 0 1 0 0 0 0 0 0
30 8971 5382 0 0 0 0 0 1 0 0 0 0 0
31 10063 4459 0 0 0 0 0 0 1 0 0 0 0
32 9793 6398 0 0 0 0 0 0 0 1 0 0 0
33 9454 4596 0 0 0 0 0 0 0 0 1 0 0
34 9759 3024 0 0 0 0 0 0 0 0 0 1 0
35 8820 1887 0 0 0 0 0 0 0 0 0 0 1
36 9403 2070 0 0 0 0 0 0 0 0 0 0 0
37 9676 1351 1 0 0 0 0 0 0 0 0 0 0
38 8642 2218 0 1 0 0 0 0 0 0 0 0 0
39 9402 2461 0 0 1 0 0 0 0 0 0 0 0
40 9610 3028 0 0 0 1 0 0 0 0 0 0 0
41 9294 4784 0 0 0 0 1 0 0 0 0 0 0
42 9448 4975 0 0 0 0 0 1 0 0 0 0 0
43 10319 4607 0 0 0 0 0 0 1 0 0 0 0
44 9548 6249 0 0 0 0 0 0 0 1 0 0 0
45 9801 4809 0 0 0 0 0 0 0 0 1 0 0
46 9596 3157 0 0 0 0 0 0 0 0 0 1 0
47 8923 1910 0 0 0 0 0 0 0 0 0 0 1
48 9746 2228 0 0 0 0 0 0 0 0 0 0 0
49 9829 1594 1 0 0 0 0 0 0 0 0 0 0
50 9125 2467 0 1 0 0 0 0 0 0 0 0 0
51 9782 2222 0 0 1 0 0 0 0 0 0 0 0
52 9441 3607 0 0 0 1 0 0 0 0 0 0 0
53 9162 4685 0 0 0 0 1 0 0 0 0 0 0
54 9915 4962 0 0 0 0 0 1 0 0 0 0 0
55 10444 5770 0 0 0 0 0 0 1 0 0 0 0
56 10209 5480 0 0 0 0 0 0 0 1 0 0 0
57 9985 5000 0 0 0 0 0 0 0 0 1 0 0
58 9842 3228 0 0 0 0 0 0 0 0 0 1 0
59 9429 1993 0 0 0 0 0 0 0 0 0 0 1
60 10132 2288 0 0 0 0 0 0 0 0 0 0 0
61 9849 1580 1 0 0 0 0 0 0 0 0 0 0
62 9172 2111 0 1 0 0 0 0 0 0 0 0 0
63 10313 2192 0 0 1 0 0 0 0 0 0 0 0
64 9819 3601 0 0 0 1 0 0 0 0 0 0 0
65 9955 4665 0 0 0 0 1 0 0 0 0 0 0
66 10048 4876 0 0 0 0 0 1 0 0 0 0 0
67 10082 5813 0 0 0 0 0 0 1 0 0 0 0
68 10541 5589 0 0 0 0 0 0 0 1 0 0 0
69 10208 5331 0 0 0 0 0 0 0 0 1 0 0
70 10233 3075 0 0 0 0 0 0 0 0 0 1 0
71 9439 2002 0 0 0 0 0 0 0 0 0 0 1
72 9963 2306 0 0 0 0 0 0 0 0 0 0 0
73 10158 1507 1 0 0 0 0 0 0 0 0 0 0
74 9225 1992 0 1 0 0 0 0 0 0 0 0 0
75 10474 2487 0 0 1 0 0 0 0 0 0 0 0
76 9757 3490 0 0 0 1 0 0 0 0 0 0 0
77 10490 4647 0 0 0 0 1 0 0 0 0 0 0
78 10281 5594 0 0 0 0 0 1 0 0 0 0 0
79 10444 5611 0 0 0 0 0 0 1 0 0 0 0
80 10640 5788 0 0 0 0 0 0 0 1 0 0 0
81 10695 6204 0 0 0 0 0 0 0 0 1 0 0
82 10786 3013 0 0 0 0 0 0 0 0 0 1 0
83 9832 1931 0 0 0 0 0 0 0 0 0 0 1
84 9747 2549 0 0 0 0 0 0 0 0 0 0 0
85 10411 1504 1 0 0 0 0 0 0 0 0 0 0
86 9511 2090 0 1 0 0 0 0 0 0 0 0 0
87 10402 2702 0 0 1 0 0 0 0 0 0 0 0
88 9701 2939 0 0 0 1 0 0 0 0 0 0 0
89 10540 4500 0 0 0 0 1 0 0 0 0 0 0
90 10112 6208 0 0 0 0 0 1 0 0 0 0 0
91 10915 6415 0 0 0 0 0 0 1 0 0 0 0
92 11183 5657 0 0 0 0 0 0 0 1 0 0 0
93 10384 5964 0 0 0 0 0 0 0 0 1 0 0
94 10834 3163 0 0 0 0 0 0 0 0 0 1 0
95 9886 1997 0 0 0 0 0 0 0 0 0 0 1
96 10216 2422 0 0 0 0 0 0 0 0 0 0 0
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) huwelijken M1 M2 M3 M4
9408.2509 0.1380 298.2272 -632.2415 245.7866 -308.7456
M5 M6 M7 M8 M9 M10
-150.9935 -429.4379 142.3794 52.8310 -237.8329 271.3407
M11
-320.0114
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-775.16 -266.33 -10.49 256.56 941.51
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 9408.2509 307.5101 30.595 < 2e-16 ***
huwelijken 0.1380 0.1171 1.178 0.24208
M1 298.2272 223.9724 1.332 0.18666
M2 -632.2415 201.3034 -3.141 0.00234 **
M3 245.7866 200.5447 1.226 0.22382
M4 -308.7456 226.7038 -1.362 0.17692
M5 -150.9935 333.5104 -0.453 0.65192
M6 -429.4379 406.8720 -1.055 0.29428
M7 142.3794 414.2315 0.344 0.73193
M8 52.8310 456.3590 0.116 0.90812
M9 -237.8329 418.7617 -0.568 0.57161
M10 271.3407 217.3278 1.249 0.21535
M11 -320.0114 206.4267 -1.550 0.12489
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 400.5 on 83 degrees of freedom
Multiple R-squared: 0.4673, Adjusted R-squared: 0.3903
F-statistic: 6.067 on 12 and 83 DF, p-value: 1.658e-07
> 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.2637129232 0.5274258463 0.73628708
[2,] 0.1560139541 0.3120279083 0.84398605
[3,] 0.0807598001 0.1615196003 0.91924020
[4,] 0.0366445899 0.0732891797 0.96335541
[5,] 0.0199034936 0.0398069873 0.98009651
[6,] 0.0150799651 0.0301599301 0.98492003
[7,] 0.0064888993 0.0129777986 0.99351110
[8,] 0.0029260678 0.0058521357 0.99707393
[9,] 0.0015287957 0.0030575915 0.99847120
[10,] 0.0008627938 0.0017255875 0.99913721
[11,] 0.0003650518 0.0007301036 0.99963495
[12,] 0.0001687943 0.0003375887 0.99983121
[13,] 0.0002771272 0.0005542543 0.99972287
[14,] 0.0024157708 0.0048315417 0.99758423
[15,] 0.0060981004 0.0121962008 0.99390190
[16,] 0.0033055682 0.0066111364 0.99669443
[17,] 0.0040171174 0.0080342349 0.99598288
[18,] 0.0025995467 0.0051990933 0.99740045
[19,] 0.0016711501 0.0033423002 0.99832885
[20,] 0.0024923694 0.0049847387 0.99750763
[21,] 0.0016596440 0.0033192880 0.99834036
[22,] 0.0010088463 0.0020176926 0.99899115
[23,] 0.0012067758 0.0024135516 0.99879322
[24,] 0.0028206150 0.0056412299 0.99717939
[25,] 0.0020676947 0.0041353894 0.99793231
[26,] 0.0073721991 0.0147443983 0.99262780
[27,] 0.0062980337 0.0125960675 0.99370197
[28,] 0.0045060334 0.0090120669 0.99549397
[29,] 0.0235838818 0.0471677635 0.97641612
[30,] 0.0214141738 0.0428283477 0.97858583
[31,] 0.0352692647 0.0705385295 0.96473074
[32,] 0.0464468475 0.0928936949 0.95355315
[33,] 0.0478840688 0.0957681375 0.95211593
[34,] 0.0383685806 0.0767371611 0.96163142
[35,] 0.0341552253 0.0683104507 0.96584477
[36,] 0.0350373991 0.0700747982 0.96496260
[37,] 0.0297454988 0.0594909976 0.97025450
[38,] 0.2595448477 0.5190896954 0.74045515
[39,] 0.2869124917 0.5738249834 0.71308751
[40,] 0.2710389095 0.5420778190 0.72896109
[41,] 0.3236268094 0.6472536188 0.67637319
[42,] 0.3336877430 0.6673754860 0.66631226
[43,] 0.5568655635 0.8862688730 0.44313444
[44,] 0.5727712325 0.8544575350 0.42722877
[45,] 0.6376788675 0.7246422650 0.36232113
[46,] 0.6644297470 0.6711405060 0.33557025
[47,] 0.6280140048 0.7439719903 0.37198600
[48,] 0.6215635913 0.7568728173 0.37843641
[49,] 0.5877626460 0.8244747081 0.41223735
[50,] 0.7048801865 0.5902396271 0.29511981
[51,] 0.6975194450 0.6049611100 0.30248055
[52,] 0.7742056840 0.4515886319 0.22579432
[53,] 0.8031877485 0.3936245030 0.19681225
[54,] 0.7927110538 0.4145778923 0.20728895
[55,] 0.9014266626 0.1971466748 0.09857334
[56,] 0.9332353937 0.1335292126 0.06676461
[57,] 0.9011845949 0.1976308102 0.09881541
[58,] 0.8818746759 0.2362506482 0.11812532
[59,] 0.8560475013 0.2879049974 0.14395250
[60,] 0.8156195655 0.3687608691 0.18438043
[61,] 0.7318192386 0.5363615228 0.26818076
[62,] 0.6709510728 0.6580978545 0.32904893
[63,] 0.6665125492 0.6669749016 0.33348745
[64,] 0.5947426776 0.8105146448 0.40525732
[65,] 0.6982910154 0.6034179692 0.30170898
> postscript(file="/var/www/rcomp/tmp/1u4h11292150293.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/rcomp/tmp/2nvy41292150293.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/rcomp/tmp/3nvy41292150293.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/rcomp/tmp/4nvy41292150293.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/rcomp/tmp/5nvy41292150293.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
-155.308508 248.939859 -54.681880 -273.247844 271.283551 -223.072329
7 8 9 10 11 12
-402.847011 -42.700843 -469.496806 -192.005408 -23.418587 -342.685998
13 14 15 16 17 18
13.072410 -374.058084 -174.548039 131.978628 209.328152 -264.913769
19 20 21 22 23 24
-52.702237 -254.183341 -633.688997 -206.210626 -159.315635 -502.342825
25 26 27 28 29 30
-208.747091 -169.163778 -209.297534 -390.189530 -423.601842 -750.285089
31 32 33 34 35 36
-102.770186 -550.715926 -350.457726 -337.766559 -528.559970 -290.817097
37 38 39 40 41 42
-216.854842 -439.992878 -591.543925 92.767925 -623.232580 -217.137535
43 44 45 46 47 48
132.812522 -775.160679 -32.842072 -519.114531 -428.732927 30.386065
49 50 51 52 53 54
-97.377828 8.656408 -178.572757 -156.107832 -741.575067 251.655875
55 56 57 58 59 60
97.371234 -8.073533 124.808584 -282.909313 65.816835 408.108784
61 62 63 64 65 66
-75.446463 104.768273 356.565883 222.719896 54.184027 396.519977
67 68 69 70 71 72
-270.560817 308.889407 302.145587 129.197753 74.575243 236.625600
73 74 75 76 77 78
243.624229 174.184879 476.869253 176.032865 591.667211 530.468520
79 80 81 82 83 84
119.306027 380.436426 668.711154 690.750943 477.370025 -12.897386
85 86 87 88 89 90
497.038093 446.665321 375.208998 196.045891 661.946548 276.764349
91 92 93 94 95 96
479.390468 941.508489 390.820277 718.057741 522.265016 473.622857
> postscript(file="/var/www/rcomp/tmp/6gmxp1292150293.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 -155.308508 NA
1 248.939859 -155.308508
2 -54.681880 248.939859
3 -273.247844 -54.681880
4 271.283551 -273.247844
5 -223.072329 271.283551
6 -402.847011 -223.072329
7 -42.700843 -402.847011
8 -469.496806 -42.700843
9 -192.005408 -469.496806
10 -23.418587 -192.005408
11 -342.685998 -23.418587
12 13.072410 -342.685998
13 -374.058084 13.072410
14 -174.548039 -374.058084
15 131.978628 -174.548039
16 209.328152 131.978628
17 -264.913769 209.328152
18 -52.702237 -264.913769
19 -254.183341 -52.702237
20 -633.688997 -254.183341
21 -206.210626 -633.688997
22 -159.315635 -206.210626
23 -502.342825 -159.315635
24 -208.747091 -502.342825
25 -169.163778 -208.747091
26 -209.297534 -169.163778
27 -390.189530 -209.297534
28 -423.601842 -390.189530
29 -750.285089 -423.601842
30 -102.770186 -750.285089
31 -550.715926 -102.770186
32 -350.457726 -550.715926
33 -337.766559 -350.457726
34 -528.559970 -337.766559
35 -290.817097 -528.559970
36 -216.854842 -290.817097
37 -439.992878 -216.854842
38 -591.543925 -439.992878
39 92.767925 -591.543925
40 -623.232580 92.767925
41 -217.137535 -623.232580
42 132.812522 -217.137535
43 -775.160679 132.812522
44 -32.842072 -775.160679
45 -519.114531 -32.842072
46 -428.732927 -519.114531
47 30.386065 -428.732927
48 -97.377828 30.386065
49 8.656408 -97.377828
50 -178.572757 8.656408
51 -156.107832 -178.572757
52 -741.575067 -156.107832
53 251.655875 -741.575067
54 97.371234 251.655875
55 -8.073533 97.371234
56 124.808584 -8.073533
57 -282.909313 124.808584
58 65.816835 -282.909313
59 408.108784 65.816835
60 -75.446463 408.108784
61 104.768273 -75.446463
62 356.565883 104.768273
63 222.719896 356.565883
64 54.184027 222.719896
65 396.519977 54.184027
66 -270.560817 396.519977
67 308.889407 -270.560817
68 302.145587 308.889407
69 129.197753 302.145587
70 74.575243 129.197753
71 236.625600 74.575243
72 243.624229 236.625600
73 174.184879 243.624229
74 476.869253 174.184879
75 176.032865 476.869253
76 591.667211 176.032865
77 530.468520 591.667211
78 119.306027 530.468520
79 380.436426 119.306027
80 668.711154 380.436426
81 690.750943 668.711154
82 477.370025 690.750943
83 -12.897386 477.370025
84 497.038093 -12.897386
85 446.665321 497.038093
86 375.208998 446.665321
87 196.045891 375.208998
88 661.946548 196.045891
89 276.764349 661.946548
90 479.390468 276.764349
91 941.508489 479.390468
92 390.820277 941.508489
93 718.057741 390.820277
94 522.265016 718.057741
95 473.622857 522.265016
96 NA 473.622857
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 248.939859 -155.308508
[2,] -54.681880 248.939859
[3,] -273.247844 -54.681880
[4,] 271.283551 -273.247844
[5,] -223.072329 271.283551
[6,] -402.847011 -223.072329
[7,] -42.700843 -402.847011
[8,] -469.496806 -42.700843
[9,] -192.005408 -469.496806
[10,] -23.418587 -192.005408
[11,] -342.685998 -23.418587
[12,] 13.072410 -342.685998
[13,] -374.058084 13.072410
[14,] -174.548039 -374.058084
[15,] 131.978628 -174.548039
[16,] 209.328152 131.978628
[17,] -264.913769 209.328152
[18,] -52.702237 -264.913769
[19,] -254.183341 -52.702237
[20,] -633.688997 -254.183341
[21,] -206.210626 -633.688997
[22,] -159.315635 -206.210626
[23,] -502.342825 -159.315635
[24,] -208.747091 -502.342825
[25,] -169.163778 -208.747091
[26,] -209.297534 -169.163778
[27,] -390.189530 -209.297534
[28,] -423.601842 -390.189530
[29,] -750.285089 -423.601842
[30,] -102.770186 -750.285089
[31,] -550.715926 -102.770186
[32,] -350.457726 -550.715926
[33,] -337.766559 -350.457726
[34,] -528.559970 -337.766559
[35,] -290.817097 -528.559970
[36,] -216.854842 -290.817097
[37,] -439.992878 -216.854842
[38,] -591.543925 -439.992878
[39,] 92.767925 -591.543925
[40,] -623.232580 92.767925
[41,] -217.137535 -623.232580
[42,] 132.812522 -217.137535
[43,] -775.160679 132.812522
[44,] -32.842072 -775.160679
[45,] -519.114531 -32.842072
[46,] -428.732927 -519.114531
[47,] 30.386065 -428.732927
[48,] -97.377828 30.386065
[49,] 8.656408 -97.377828
[50,] -178.572757 8.656408
[51,] -156.107832 -178.572757
[52,] -741.575067 -156.107832
[53,] 251.655875 -741.575067
[54,] 97.371234 251.655875
[55,] -8.073533 97.371234
[56,] 124.808584 -8.073533
[57,] -282.909313 124.808584
[58,] 65.816835 -282.909313
[59,] 408.108784 65.816835
[60,] -75.446463 408.108784
[61,] 104.768273 -75.446463
[62,] 356.565883 104.768273
[63,] 222.719896 356.565883
[64,] 54.184027 222.719896
[65,] 396.519977 54.184027
[66,] -270.560817 396.519977
[67,] 308.889407 -270.560817
[68,] 302.145587 308.889407
[69,] 129.197753 302.145587
[70,] 74.575243 129.197753
[71,] 236.625600 74.575243
[72,] 243.624229 236.625600
[73,] 174.184879 243.624229
[74,] 476.869253 174.184879
[75,] 176.032865 476.869253
[76,] 591.667211 176.032865
[77,] 530.468520 591.667211
[78,] 119.306027 530.468520
[79,] 380.436426 119.306027
[80,] 668.711154 380.436426
[81,] 690.750943 668.711154
[82,] 477.370025 690.750943
[83,] -12.897386 477.370025
[84,] 497.038093 -12.897386
[85,] 446.665321 497.038093
[86,] 375.208998 446.665321
[87,] 196.045891 375.208998
[88,] 661.946548 196.045891
[89,] 276.764349 661.946548
[90,] 479.390468 276.764349
[91,] 941.508489 479.390468
[92,] 390.820277 941.508489
[93,] 718.057741 390.820277
[94,] 522.265016 718.057741
[95,] 473.622857 522.265016
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 248.939859 -155.308508
2 -54.681880 248.939859
3 -273.247844 -54.681880
4 271.283551 -273.247844
5 -223.072329 271.283551
6 -402.847011 -223.072329
7 -42.700843 -402.847011
8 -469.496806 -42.700843
9 -192.005408 -469.496806
10 -23.418587 -192.005408
11 -342.685998 -23.418587
12 13.072410 -342.685998
13 -374.058084 13.072410
14 -174.548039 -374.058084
15 131.978628 -174.548039
16 209.328152 131.978628
17 -264.913769 209.328152
18 -52.702237 -264.913769
19 -254.183341 -52.702237
20 -633.688997 -254.183341
21 -206.210626 -633.688997
22 -159.315635 -206.210626
23 -502.342825 -159.315635
24 -208.747091 -502.342825
25 -169.163778 -208.747091
26 -209.297534 -169.163778
27 -390.189530 -209.297534
28 -423.601842 -390.189530
29 -750.285089 -423.601842
30 -102.770186 -750.285089
31 -550.715926 -102.770186
32 -350.457726 -550.715926
33 -337.766559 -350.457726
34 -528.559970 -337.766559
35 -290.817097 -528.559970
36 -216.854842 -290.817097
37 -439.992878 -216.854842
38 -591.543925 -439.992878
39 92.767925 -591.543925
40 -623.232580 92.767925
41 -217.137535 -623.232580
42 132.812522 -217.137535
43 -775.160679 132.812522
44 -32.842072 -775.160679
45 -519.114531 -32.842072
46 -428.732927 -519.114531
47 30.386065 -428.732927
48 -97.377828 30.386065
49 8.656408 -97.377828
50 -178.572757 8.656408
51 -156.107832 -178.572757
52 -741.575067 -156.107832
53 251.655875 -741.575067
54 97.371234 251.655875
55 -8.073533 97.371234
56 124.808584 -8.073533
57 -282.909313 124.808584
58 65.816835 -282.909313
59 408.108784 65.816835
60 -75.446463 408.108784
61 104.768273 -75.446463
62 356.565883 104.768273
63 222.719896 356.565883
64 54.184027 222.719896
65 396.519977 54.184027
66 -270.560817 396.519977
67 308.889407 -270.560817
68 302.145587 308.889407
69 129.197753 302.145587
70 74.575243 129.197753
71 236.625600 74.575243
72 243.624229 236.625600
73 174.184879 243.624229
74 476.869253 174.184879
75 176.032865 476.869253
76 591.667211 176.032865
77 530.468520 591.667211
78 119.306027 530.468520
79 380.436426 119.306027
80 668.711154 380.436426
81 690.750943 668.711154
82 477.370025 690.750943
83 -12.897386 477.370025
84 497.038093 -12.897386
85 446.665321 497.038093
86 375.208998 446.665321
87 196.045891 375.208998
88 661.946548 196.045891
89 276.764349 661.946548
90 479.390468 276.764349
91 941.508489 479.390468
92 390.820277 941.508489
93 718.057741 390.820277
94 522.265016 718.057741
95 473.622857 522.265016
> 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/rcomp/tmp/7rvfs1292150293.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/rcomp/tmp/8rvfs1292150293.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/rcomp/tmp/9rvfs1292150293.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/rcomp/tmp/1015ev1292150293.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/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/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/rcomp/tmp/11n5cj1292150293.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/rcomp/tmp/12q6b71292150293.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/rcomp/tmp/13mf8x1292150293.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/rcomp/tmp/14pg7l1292150293.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/rcomp/tmp/15byn91292150293.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/rcomp/tmp/167ql01292150293.tab")
+ }
>
> try(system("convert tmp/1u4h11292150293.ps tmp/1u4h11292150293.png",intern=TRUE))
character(0)
> try(system("convert tmp/2nvy41292150293.ps tmp/2nvy41292150293.png",intern=TRUE))
character(0)
> try(system("convert tmp/3nvy41292150293.ps tmp/3nvy41292150293.png",intern=TRUE))
character(0)
> try(system("convert tmp/4nvy41292150293.ps tmp/4nvy41292150293.png",intern=TRUE))
character(0)
> try(system("convert tmp/5nvy41292150293.ps tmp/5nvy41292150293.png",intern=TRUE))
character(0)
> try(system("convert tmp/6gmxp1292150293.ps tmp/6gmxp1292150293.png",intern=TRUE))
character(0)
> try(system("convert tmp/7rvfs1292150293.ps tmp/7rvfs1292150293.png",intern=TRUE))
character(0)
> try(system("convert tmp/8rvfs1292150293.ps tmp/8rvfs1292150293.png",intern=TRUE))
character(0)
> try(system("convert tmp/9rvfs1292150293.ps tmp/9rvfs1292150293.png",intern=TRUE))
character(0)
> try(system("convert tmp/1015ev1292150293.ps tmp/1015ev1292150293.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.560 1.690 5.215