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(5393
+ ,552486
+ ,3.90
+ ,3.0
+ ,628232
+ ,5147
+ ,516610
+ ,3.90
+ ,2.2
+ ,612117
+ ,4846
+ ,487587
+ ,3.88
+ ,2.3
+ ,595404
+ ,3995
+ ,403620
+ ,3.89
+ ,2.8
+ ,597141
+ ,4491
+ ,459427
+ ,3.89
+ ,2.8
+ ,593408
+ ,4676
+ ,473058
+ ,3.93
+ ,2.8
+ ,590072
+ ,5461
+ ,583054
+ ,3.94
+ ,2.2
+ ,579799
+ ,4758
+ ,509448
+ ,3.97
+ ,2.6
+ ,574205
+ ,5302
+ ,551582
+ ,4.00
+ ,2.8
+ ,572775
+ ,5066
+ ,524752
+ ,4.04
+ ,2.5
+ ,572942
+ ,3491
+ ,370725
+ ,4.18
+ ,2.4
+ ,619567
+ ,4944
+ ,531443
+ ,4.32
+ ,2.3
+ ,625809
+ ,5148
+ ,537833
+ ,4.37
+ ,1.9
+ ,619916
+ ,5351
+ ,551410
+ ,4.40
+ ,1.7
+ ,587625
+ ,5178
+ ,520983
+ ,4.38
+ ,2.0
+ ,565742
+ ,4025
+ ,395542
+ ,4.36
+ ,2.1
+ ,557274
+ ,4449
+ ,442878
+ ,4.36
+ ,1.7
+ ,560576
+ ,4594
+ ,454919
+ ,4.40
+ ,1.8
+ ,548854
+ ,4603
+ ,488905
+ ,4.41
+ ,1.8
+ ,531673
+ ,4911
+ ,496085
+ ,4.43
+ ,1.8
+ ,525919
+ ,5236
+ ,540146
+ ,4.42
+ ,1.3
+ ,511038
+ ,4652
+ ,496529
+ ,4.46
+ ,1.3
+ ,498662
+ ,3479
+ ,372656
+ ,4.61
+ ,1.3
+ ,555362
+ ,4556
+ ,486704
+ ,4.78
+ ,1.2
+ ,564591
+ ,4815
+ ,495334
+ ,4.88
+ ,1.4
+ ,541657
+ ,4949
+ ,504697
+ ,4.95
+ ,2.2
+ ,527070
+ ,4499
+ ,464856
+ ,4.95
+ ,2.9
+ ,509846
+ ,3865
+ ,388472
+ ,4.93
+ ,3.1
+ ,514258
+ ,3657
+ ,377508
+ ,4.93
+ ,3.5
+ ,516922
+ ,4814
+ ,468895
+ ,4.91
+ ,3.6
+ ,507561
+ ,4614
+ ,471295
+ ,4.88
+ ,4.4
+ ,492622
+ ,4539
+ ,482956
+ ,4.83
+ ,4.1
+ ,490243
+ ,4492
+ ,483404
+ ,4.83
+ ,5.1
+ ,469357
+ ,4779
+ ,495548
+ ,4.85
+ ,5.8
+ ,477580
+ ,3193
+ ,333806
+ ,4.99
+ ,5.9
+ ,528379
+ ,3894
+ ,411611
+ ,5.14
+ ,5.4
+ ,533590
+ ,4531
+ ,496215
+ ,5.26
+ ,5.5
+ ,517945
+ ,4008
+ ,433542
+ ,5.33
+ ,4.8
+ ,506174
+ ,3764
+ ,409819
+ ,5.28
+ ,3.2
+ ,501866
+ ,3290
+ ,339270
+ ,4.99
+ ,2.7
+ ,516141
+ ,3644
+ ,365092
+ ,4.75
+ ,2.1
+ ,528222
+ ,3438
+ ,387851
+ ,4.63
+ ,1.9
+ ,532638
+ ,3833
+ ,408171
+ ,4.52
+ ,0.6
+ ,536322
+ ,3922
+ ,427587
+ ,4.50
+ ,0.7
+ ,536535
+ ,3524
+ ,377805
+ ,4.48
+ ,-0.2
+ ,523597
+ ,3493
+ ,376222
+ ,4.49
+ ,-1.0
+ ,536214
+ ,2814
+ ,300606
+ ,4.57
+ ,-1.7
+ ,586570
+ ,3899
+ ,424611
+ ,4.64
+ ,-0.7
+ ,596594
+ ,3653
+ ,404393
+ ,4.62
+ ,-1.0
+ ,580523
+ ,3969
+ ,422701
+ ,4.55
+ ,-0.9
+ ,564478
+ ,3427
+ ,369704
+ ,4.47
+ ,0.0
+ ,557560
+ ,3067
+ ,320685
+ ,4.43
+ ,0.3
+ ,575093
+ ,3301
+ ,344674
+ ,4.45
+ ,0.8
+ ,580112
+ ,3211
+ ,319302
+ ,4.41
+ ,0.8
+ ,574761
+ ,3382
+ ,368391
+ ,4.32
+ ,1.9
+ ,563250
+ ,3613
+ ,395375
+ ,4.24
+ ,2.1
+ ,551531
+ ,3783
+ ,420926
+ ,4.16
+ ,2.5
+ ,537034
+ ,3971
+ ,434358
+ ,4.03
+ ,2.7
+ ,544686
+ ,2842
+ ,315828
+ ,4.01
+ ,2.4
+ ,600991
+ ,4161
+ ,451722
+ ,3.98
+ ,2.4
+ ,604378)
+ ,dim=c(5
+ ,60)
+ ,dimnames=list(c('Nieuwe_woningen'
+ ,'Bewoonbare_opp'
+ ,'Rentevoet'
+ ,'Inflatie'
+ ,'Werkloosheid')
+ ,1:60))
> y <- array(NA,dim=c(5,60),dimnames=list(c('Nieuwe_woningen','Bewoonbare_opp','Rentevoet','Inflatie','Werkloosheid'),1:60))
> 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
Nieuwe_woningen Bewoonbare_opp Rentevoet Inflatie Werkloosheid M1 M2 M3 M4
1 5393 552486 3.90 3.0 628232 1 0 0 0
2 5147 516610 3.90 2.2 612117 0 1 0 0
3 4846 487587 3.88 2.3 595404 0 0 1 0
4 3995 403620 3.89 2.8 597141 0 0 0 1
5 4491 459427 3.89 2.8 593408 0 0 0 0
6 4676 473058 3.93 2.8 590072 0 0 0 0
7 5461 583054 3.94 2.2 579799 0 0 0 0
8 4758 509448 3.97 2.6 574205 0 0 0 0
9 5302 551582 4.00 2.8 572775 0 0 0 0
10 5066 524752 4.04 2.5 572942 0 0 0 0
11 3491 370725 4.18 2.4 619567 0 0 0 0
12 4944 531443 4.32 2.3 625809 0 0 0 0
13 5148 537833 4.37 1.9 619916 1 0 0 0
14 5351 551410 4.40 1.7 587625 0 1 0 0
15 5178 520983 4.38 2.0 565742 0 0 1 0
16 4025 395542 4.36 2.1 557274 0 0 0 1
17 4449 442878 4.36 1.7 560576 0 0 0 0
18 4594 454919 4.40 1.8 548854 0 0 0 0
19 4603 488905 4.41 1.8 531673 0 0 0 0
20 4911 496085 4.43 1.8 525919 0 0 0 0
21 5236 540146 4.42 1.3 511038 0 0 0 0
22 4652 496529 4.46 1.3 498662 0 0 0 0
23 3479 372656 4.61 1.3 555362 0 0 0 0
24 4556 486704 4.78 1.2 564591 0 0 0 0
25 4815 495334 4.88 1.4 541657 1 0 0 0
26 4949 504697 4.95 2.2 527070 0 1 0 0
27 4499 464856 4.95 2.9 509846 0 0 1 0
28 3865 388472 4.93 3.1 514258 0 0 0 1
29 3657 377508 4.93 3.5 516922 0 0 0 0
30 4814 468895 4.91 3.6 507561 0 0 0 0
31 4614 471295 4.88 4.4 492622 0 0 0 0
32 4539 482956 4.83 4.1 490243 0 0 0 0
33 4492 483404 4.83 5.1 469357 0 0 0 0
34 4779 495548 4.85 5.8 477580 0 0 0 0
35 3193 333806 4.99 5.9 528379 0 0 0 0
36 3894 411611 5.14 5.4 533590 0 0 0 0
37 4531 496215 5.26 5.5 517945 1 0 0 0
38 4008 433542 5.33 4.8 506174 0 1 0 0
39 3764 409819 5.28 3.2 501866 0 0 1 0
40 3290 339270 4.99 2.7 516141 0 0 0 1
41 3644 365092 4.75 2.1 528222 0 0 0 0
42 3438 387851 4.63 1.9 532638 0 0 0 0
43 3833 408171 4.52 0.6 536322 0 0 0 0
44 3922 427587 4.50 0.7 536535 0 0 0 0
45 3524 377805 4.48 -0.2 523597 0 0 0 0
46 3493 376222 4.49 -1.0 536214 0 0 0 0
47 2814 300606 4.57 -1.7 586570 0 0 0 0
48 3899 424611 4.64 -0.7 596594 0 0 0 0
49 3653 404393 4.62 -1.0 580523 1 0 0 0
50 3969 422701 4.55 -0.9 564478 0 1 0 0
51 3427 369704 4.47 0.0 557560 0 0 1 0
52 3067 320685 4.43 0.3 575093 0 0 0 1
53 3301 344674 4.45 0.8 580112 0 0 0 0
54 3211 319302 4.41 0.8 574761 0 0 0 0
55 3382 368391 4.32 1.9 563250 0 0 0 0
56 3613 395375 4.24 2.1 551531 0 0 0 0
57 3783 420926 4.16 2.5 537034 0 0 0 0
58 3971 434358 4.03 2.7 544686 0 0 0 0
59 2842 315828 4.01 2.4 600991 0 0 0 0
60 4161 451722 3.98 2.4 604378 0 0 0 0
M5 M6 M7 M8 M9 M10 M11 t
1 0 0 0 0 0 0 0 1
2 0 0 0 0 0 0 0 2
3 0 0 0 0 0 0 0 3
4 0 0 0 0 0 0 0 4
5 1 0 0 0 0 0 0 5
6 0 1 0 0 0 0 0 6
7 0 0 1 0 0 0 0 7
8 0 0 0 1 0 0 0 8
9 0 0 0 0 1 0 0 9
10 0 0 0 0 0 1 0 10
11 0 0 0 0 0 0 1 11
12 0 0 0 0 0 0 0 12
13 0 0 0 0 0 0 0 13
14 0 0 0 0 0 0 0 14
15 0 0 0 0 0 0 0 15
16 0 0 0 0 0 0 0 16
17 1 0 0 0 0 0 0 17
18 0 1 0 0 0 0 0 18
19 0 0 1 0 0 0 0 19
20 0 0 0 1 0 0 0 20
21 0 0 0 0 1 0 0 21
22 0 0 0 0 0 1 0 22
23 0 0 0 0 0 0 1 23
24 0 0 0 0 0 0 0 24
25 0 0 0 0 0 0 0 25
26 0 0 0 0 0 0 0 26
27 0 0 0 0 0 0 0 27
28 0 0 0 0 0 0 0 28
29 1 0 0 0 0 0 0 29
30 0 1 0 0 0 0 0 30
31 0 0 1 0 0 0 0 31
32 0 0 0 1 0 0 0 32
33 0 0 0 0 1 0 0 33
34 0 0 0 0 0 1 0 34
35 0 0 0 0 0 0 1 35
36 0 0 0 0 0 0 0 36
37 0 0 0 0 0 0 0 37
38 0 0 0 0 0 0 0 38
39 0 0 0 0 0 0 0 39
40 0 0 0 0 0 0 0 40
41 1 0 0 0 0 0 0 41
42 0 1 0 0 0 0 0 42
43 0 0 1 0 0 0 0 43
44 0 0 0 1 0 0 0 44
45 0 0 0 0 1 0 0 45
46 0 0 0 0 0 1 0 46
47 0 0 0 0 0 0 1 47
48 0 0 0 0 0 0 0 48
49 0 0 0 0 0 0 0 49
50 0 0 0 0 0 0 0 50
51 0 0 0 0 0 0 0 51
52 0 0 0 0 0 0 0 52
53 1 0 0 0 0 0 0 53
54 0 1 0 0 0 0 0 54
55 0 0 1 0 0 0 0 55
56 0 0 0 1 0 0 0 56
57 0 0 0 0 1 0 0 57
58 0 0 0 0 0 1 0 58
59 0 0 0 0 0 0 1 59
60 0 0 0 0 0 0 0 60
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Bewoonbare_opp Rentevoet Inflatie Werkloosheid
844.662453 0.009800 -59.355052 -6.879358 -0.001117
M1 M2 M3 M4 M5
17.023684 89.681963 79.962178 186.938946 173.949223
M6 M7 M8 M9 M10
184.728747 -16.184774 -32.320235 -47.213120 -24.499993
M11 t
55.729790 -3.740587
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-352.62 -59.18 12.85 63.80 184.59
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 8.447e+02 1.440e+03 0.587 0.5605
Bewoonbare_opp 9.800e-03 6.014e-04 16.296 <2e-16 ***
Rentevoet -5.936e+01 1.173e+02 -0.506 0.6154
Inflatie -6.879e+00 1.099e+01 -0.626 0.5345
Werkloosheid -1.117e-03 1.364e-03 -0.819 0.4175
M1 1.702e+01 7.039e+01 0.242 0.8101
M2 8.968e+01 7.616e+01 1.178 0.2455
M3 7.996e+01 9.165e+01 0.873 0.3878
M4 1.869e+02 1.169e+02 1.600 0.1170
M5 1.739e+02 1.048e+02 1.661 0.1041
M6 1.847e+02 1.023e+02 1.806 0.0780 .
M7 -1.618e+01 1.053e+02 -0.154 0.8786
M8 -3.232e+01 1.119e+02 -0.289 0.7742
M9 -4.721e+01 1.255e+02 -0.376 0.7086
M10 -2.450e+01 1.231e+02 -0.199 0.8432
M11 5.573e+01 1.081e+02 0.516 0.6087
t -3.741e+00 1.779e+00 -2.103 0.0413 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 109.1 on 43 degrees of freedom
Multiple R-squared: 0.9829, Adjusted R-squared: 0.9766
F-statistic: 154.7 on 16 and 43 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.06009478 0.1201896 0.9399052
[2,] 0.36646657 0.7329331 0.6335334
[3,] 0.53279866 0.9344027 0.4672013
[4,] 0.44000400 0.8800080 0.5599960
[5,] 0.43328030 0.8665606 0.5667197
[6,] 0.36048958 0.7209792 0.6395104
[7,] 0.26090644 0.5218129 0.7390936
[8,] 0.20215464 0.4043093 0.7978454
[9,] 0.13262991 0.2652598 0.8673701
[10,] 0.12574322 0.2514864 0.8742568
[11,] 0.24365733 0.4873147 0.7563427
[12,] 0.22289689 0.4457938 0.7771031
[13,] 0.35038615 0.7007723 0.6496139
[14,] 0.33709719 0.6741944 0.6629028
[15,] 0.41819133 0.8363827 0.5818087
[16,] 0.35356243 0.7071249 0.6464376
[17,] 0.25006334 0.5001267 0.7499367
[18,] 0.44788241 0.8957648 0.5521176
[19,] 0.39750654 0.7950131 0.6024935
[20,] 0.32867643 0.6573529 0.6713236
[21,] 0.21714890 0.4342978 0.7828511
> postscript(file="/var/www/html/freestat/rcomp/tmp/1dnaa1296731106.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/2eu7x1296731106.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/360sy1296731106.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/4zv381296731106.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/5kdc91296731106.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 = 60
Frequency = 1
1 2 3 4 5 6
74.282987 87.451301 65.175464 -60.218054 -98.559918 -55.532059
7 8 9 10 11 12
-158.834126 -122.338357 28.945115 37.401338 -44.956336 -92.918832
13 14 15 16 17 18
28.813001 -5.806586 109.280197 72.402299 50.178460 60.113101
19 20 21 22 23 24
-77.882102 174.393224 65.584355 -121.388672 -84.722901 -46.207373
25 26 27 28 29 30
96.641604 63.339092 2.822993 19.258863 -58.838361 184.586431
31 32 33 34 35 36
152.764061 -24.324514 -73.521417 90.679481 78.964379 88.233513
37 38 39 40 41 42
-126.819407 -118.349682 -135.190396 -25.766376 87.026042 -352.617305
43 44 45 46 47 48
36.543377 -45.117212 41.552090 16.269863 57.967266 9.424381
49 50 51 52 53 54
-72.918185 -26.634126 -42.088259 -5.676732 20.193777 163.449832
55 56 57 58 59 60
47.408790 17.386859 -62.560144 -22.962009 -7.252408 41.468311
> postscript(file="/var/www/html/freestat/rcomp/tmp/6txy21296731106.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 = 60
Frequency = 1
lag(myerror, k = 1) myerror
0 74.282987 NA
1 87.451301 74.282987
2 65.175464 87.451301
3 -60.218054 65.175464
4 -98.559918 -60.218054
5 -55.532059 -98.559918
6 -158.834126 -55.532059
7 -122.338357 -158.834126
8 28.945115 -122.338357
9 37.401338 28.945115
10 -44.956336 37.401338
11 -92.918832 -44.956336
12 28.813001 -92.918832
13 -5.806586 28.813001
14 109.280197 -5.806586
15 72.402299 109.280197
16 50.178460 72.402299
17 60.113101 50.178460
18 -77.882102 60.113101
19 174.393224 -77.882102
20 65.584355 174.393224
21 -121.388672 65.584355
22 -84.722901 -121.388672
23 -46.207373 -84.722901
24 96.641604 -46.207373
25 63.339092 96.641604
26 2.822993 63.339092
27 19.258863 2.822993
28 -58.838361 19.258863
29 184.586431 -58.838361
30 152.764061 184.586431
31 -24.324514 152.764061
32 -73.521417 -24.324514
33 90.679481 -73.521417
34 78.964379 90.679481
35 88.233513 78.964379
36 -126.819407 88.233513
37 -118.349682 -126.819407
38 -135.190396 -118.349682
39 -25.766376 -135.190396
40 87.026042 -25.766376
41 -352.617305 87.026042
42 36.543377 -352.617305
43 -45.117212 36.543377
44 41.552090 -45.117212
45 16.269863 41.552090
46 57.967266 16.269863
47 9.424381 57.967266
48 -72.918185 9.424381
49 -26.634126 -72.918185
50 -42.088259 -26.634126
51 -5.676732 -42.088259
52 20.193777 -5.676732
53 163.449832 20.193777
54 47.408790 163.449832
55 17.386859 47.408790
56 -62.560144 17.386859
57 -22.962009 -62.560144
58 -7.252408 -22.962009
59 41.468311 -7.252408
60 NA 41.468311
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 87.451301 74.282987
[2,] 65.175464 87.451301
[3,] -60.218054 65.175464
[4,] -98.559918 -60.218054
[5,] -55.532059 -98.559918
[6,] -158.834126 -55.532059
[7,] -122.338357 -158.834126
[8,] 28.945115 -122.338357
[9,] 37.401338 28.945115
[10,] -44.956336 37.401338
[11,] -92.918832 -44.956336
[12,] 28.813001 -92.918832
[13,] -5.806586 28.813001
[14,] 109.280197 -5.806586
[15,] 72.402299 109.280197
[16,] 50.178460 72.402299
[17,] 60.113101 50.178460
[18,] -77.882102 60.113101
[19,] 174.393224 -77.882102
[20,] 65.584355 174.393224
[21,] -121.388672 65.584355
[22,] -84.722901 -121.388672
[23,] -46.207373 -84.722901
[24,] 96.641604 -46.207373
[25,] 63.339092 96.641604
[26,] 2.822993 63.339092
[27,] 19.258863 2.822993
[28,] -58.838361 19.258863
[29,] 184.586431 -58.838361
[30,] 152.764061 184.586431
[31,] -24.324514 152.764061
[32,] -73.521417 -24.324514
[33,] 90.679481 -73.521417
[34,] 78.964379 90.679481
[35,] 88.233513 78.964379
[36,] -126.819407 88.233513
[37,] -118.349682 -126.819407
[38,] -135.190396 -118.349682
[39,] -25.766376 -135.190396
[40,] 87.026042 -25.766376
[41,] -352.617305 87.026042
[42,] 36.543377 -352.617305
[43,] -45.117212 36.543377
[44,] 41.552090 -45.117212
[45,] 16.269863 41.552090
[46,] 57.967266 16.269863
[47,] 9.424381 57.967266
[48,] -72.918185 9.424381
[49,] -26.634126 -72.918185
[50,] -42.088259 -26.634126
[51,] -5.676732 -42.088259
[52,] 20.193777 -5.676732
[53,] 163.449832 20.193777
[54,] 47.408790 163.449832
[55,] 17.386859 47.408790
[56,] -62.560144 17.386859
[57,] -22.962009 -62.560144
[58,] -7.252408 -22.962009
[59,] 41.468311 -7.252408
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 87.451301 74.282987
2 65.175464 87.451301
3 -60.218054 65.175464
4 -98.559918 -60.218054
5 -55.532059 -98.559918
6 -158.834126 -55.532059
7 -122.338357 -158.834126
8 28.945115 -122.338357
9 37.401338 28.945115
10 -44.956336 37.401338
11 -92.918832 -44.956336
12 28.813001 -92.918832
13 -5.806586 28.813001
14 109.280197 -5.806586
15 72.402299 109.280197
16 50.178460 72.402299
17 60.113101 50.178460
18 -77.882102 60.113101
19 174.393224 -77.882102
20 65.584355 174.393224
21 -121.388672 65.584355
22 -84.722901 -121.388672
23 -46.207373 -84.722901
24 96.641604 -46.207373
25 63.339092 96.641604
26 2.822993 63.339092
27 19.258863 2.822993
28 -58.838361 19.258863
29 184.586431 -58.838361
30 152.764061 184.586431
31 -24.324514 152.764061
32 -73.521417 -24.324514
33 90.679481 -73.521417
34 78.964379 90.679481
35 88.233513 78.964379
36 -126.819407 88.233513
37 -118.349682 -126.819407
38 -135.190396 -118.349682
39 -25.766376 -135.190396
40 87.026042 -25.766376
41 -352.617305 87.026042
42 36.543377 -352.617305
43 -45.117212 36.543377
44 41.552090 -45.117212
45 16.269863 41.552090
46 57.967266 16.269863
47 9.424381 57.967266
48 -72.918185 9.424381
49 -26.634126 -72.918185
50 -42.088259 -26.634126
51 -5.676732 -42.088259
52 20.193777 -5.676732
53 163.449832 20.193777
54 47.408790 163.449832
55 17.386859 47.408790
56 -62.560144 17.386859
57 -22.962009 -62.560144
58 -7.252408 -22.962009
59 41.468311 -7.252408
> 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/7hdrz1296731106.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/88s2f1296731106.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/96hqi1296731106.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/10gj2z1296731106.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/116xrp1296731106.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/1294db1296731106.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/13ba3q1296731106.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/142is01296731107.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/15af6u1296731107.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/16bvck1296731107.tab")
+ }
>
> try(system("convert tmp/1dnaa1296731106.ps tmp/1dnaa1296731106.png",intern=TRUE))
character(0)
> try(system("convert tmp/2eu7x1296731106.ps tmp/2eu7x1296731106.png",intern=TRUE))
character(0)
> try(system("convert tmp/360sy1296731106.ps tmp/360sy1296731106.png",intern=TRUE))
character(0)
> try(system("convert tmp/4zv381296731106.ps tmp/4zv381296731106.png",intern=TRUE))
character(0)
> try(system("convert tmp/5kdc91296731106.ps tmp/5kdc91296731106.png",intern=TRUE))
character(0)
> try(system("convert tmp/6txy21296731106.ps tmp/6txy21296731106.png",intern=TRUE))
character(0)
> try(system("convert tmp/7hdrz1296731106.ps tmp/7hdrz1296731106.png",intern=TRUE))
character(0)
> try(system("convert tmp/88s2f1296731106.ps tmp/88s2f1296731106.png",intern=TRUE))
character(0)
> try(system("convert tmp/96hqi1296731106.ps tmp/96hqi1296731106.png",intern=TRUE))
character(0)
> try(system("convert tmp/10gj2z1296731106.ps tmp/10gj2z1296731106.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.832 2.495 4.395