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 = 'No Linear Trend'
> par2 = 'Do not include Seasonal Dummies'
> par1 = '1'
> #'GNU S' R Code compiled by R2WASP v. 1.0.44 ()
> #Author: Prof. Dr. P. Wessa
> #To cite this work: AUTHOR(S), (YEAR), YOUR SOFTWARE TITLE (vNUMBER) in Free Statistics Software (v$_version), Office for Research Development and Education, URL http://www.wessa.net/rwasp_YOURPAGE.wasp/
> #Source of accompanying publication: Office for Research, Development, and Education
> #Technical description: Write here your technical program description (don't use hard returns!)
> library(lattice)
> library(lmtest)
Loading required package: zoo
Attaching package: 'zoo'
The following object(s) are masked from package:base :
as.Date.numeric
> n25 <- 25 #minimum number of obs. for Goldfeld-Quandt test
> par1 <- as.numeric(par1)
> x <- t(y)
> k <- length(x[1,])
> n <- length(x[,1])
> x1 <- cbind(x[,par1], x[,1:k!=par1])
> mycolnames <- c(colnames(x)[par1], colnames(x)[1:k!=par1])
> colnames(x1) <- mycolnames #colnames(x)[par1]
> x <- x1
> if (par3 == 'First Differences'){
+ x2 <- array(0, dim=c(n-1,k), dimnames=list(1:(n-1), paste('(1-B)',colnames(x),sep='')))
+ for (i in 1:n-1) {
+ for (j in 1:k) {
+ x2[i,j] <- x[i+1,j] - x[i,j]
+ }
+ }
+ x <- x2
+ }
> if (par2 == 'Include Monthly Dummies'){
+ x2 <- array(0, dim=c(n,11), dimnames=list(1:n, paste('M', seq(1:11), sep ='')))
+ for (i in 1:11){
+ x2[seq(i,n,12),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> if (par2 == 'Include Quarterly Dummies'){
+ x2 <- array(0, dim=c(n,3), dimnames=list(1:n, paste('Q', seq(1:3), sep ='')))
+ for (i in 1:3){
+ x2[seq(i,n,4),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> k <- length(x[1,])
> if (par3 == 'Linear Trend'){
+ x <- cbind(x, c(1:n))
+ colnames(x)[k+1] <- 't'
+ }
> x
Nieuwe_woningen Bewoonbare_opp Rentevoet Inflatie Werkloosheid
1 5393 552486 3.90 3.0 628232
2 5147 516610 3.90 2.2 612117
3 4846 487587 3.88 2.3 595404
4 3995 403620 3.89 2.8 597141
5 4491 459427 3.89 2.8 593408
6 4676 473058 3.93 2.8 590072
7 5461 583054 3.94 2.2 579799
8 4758 509448 3.97 2.6 574205
9 5302 551582 4.00 2.8 572775
10 5066 524752 4.04 2.5 572942
11 3491 370725 4.18 2.4 619567
12 4944 531443 4.32 2.3 625809
13 5148 537833 4.37 1.9 619916
14 5351 551410 4.40 1.7 587625
15 5178 520983 4.38 2.0 565742
16 4025 395542 4.36 2.1 557274
17 4449 442878 4.36 1.7 560576
18 4594 454919 4.40 1.8 548854
19 4603 488905 4.41 1.8 531673
20 4911 496085 4.43 1.8 525919
21 5236 540146 4.42 1.3 511038
22 4652 496529 4.46 1.3 498662
23 3479 372656 4.61 1.3 555362
24 4556 486704 4.78 1.2 564591
25 4815 495334 4.88 1.4 541657
26 4949 504697 4.95 2.2 527070
27 4499 464856 4.95 2.9 509846
28 3865 388472 4.93 3.1 514258
29 3657 377508 4.93 3.5 516922
30 4814 468895 4.91 3.6 507561
31 4614 471295 4.88 4.4 492622
32 4539 482956 4.83 4.1 490243
33 4492 483404 4.83 5.1 469357
34 4779 495548 4.85 5.8 477580
35 3193 333806 4.99 5.9 528379
36 3894 411611 5.14 5.4 533590
37 4531 496215 5.26 5.5 517945
38 4008 433542 5.33 4.8 506174
39 3764 409819 5.28 3.2 501866
40 3290 339270 4.99 2.7 516141
41 3644 365092 4.75 2.1 528222
42 3438 387851 4.63 1.9 532638
43 3833 408171 4.52 0.6 536322
44 3922 427587 4.50 0.7 536535
45 3524 377805 4.48 -0.2 523597
46 3493 376222 4.49 -1.0 536214
47 2814 300606 4.57 -1.7 586570
48 3899 424611 4.64 -0.7 596594
49 3653 404393 4.62 -1.0 580523
50 3969 422701 4.55 -0.9 564478
51 3427 369704 4.47 0.0 557560
52 3067 320685 4.43 0.3 575093
53 3301 344674 4.45 0.8 580112
54 3211 319302 4.41 0.8 574761
55 3382 368391 4.32 1.9 563250
56 3613 395375 4.24 2.1 551531
57 3783 420926 4.16 2.5 537034
58 3971 434358 4.03 2.7 544686
59 2842 315828 4.01 2.4 600991
60 4161 451722 3.98 2.4 604378
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Bewoonbare_opp Rentevoet Inflatie Werkloosheid
-53.684785 0.009835 -31.564587 10.488404 0.000082
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-241.95 -118.98 -19.81 95.47 331.84
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -5.368e+01 6.884e+02 -0.078 0.938
Bewoonbare_opp 9.835e-03 2.927e-04 33.597 <2e-16 ***
Rentevoet -3.156e+01 7.133e+01 -0.442 0.660
Inflatie 1.049e+01 1.300e+01 0.807 0.423
Werkloosheid 8.201e-05 7.139e-04 0.115 0.909
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 144.6 on 55 degrees of freedom
Multiple R-squared: 0.9616, Adjusted R-squared: 0.9589
F-statistic: 344.7 on 4 and 55 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.1774301 0.35486014 0.82256993
[2,] 0.2376214 0.47524288 0.76237856
[3,] 0.1295856 0.25917110 0.87041445
[4,] 0.1753269 0.35065374 0.82467313
[5,] 0.1088839 0.21776784 0.89111608
[6,] 0.1272244 0.25444882 0.87277559
[7,] 0.1840486 0.36809713 0.81595143
[8,] 0.3146196 0.62923915 0.68538043
[9,] 0.3457401 0.69148014 0.65425993
[10,] 0.3345423 0.66908453 0.66545774
[11,] 0.3712304 0.74246085 0.62876958
[12,] 0.4714035 0.94280695 0.52859653
[13,] 0.4699718 0.93994350 0.53002825
[14,] 0.4382211 0.87644223 0.56177888
[15,] 0.5176821 0.96463571 0.48231786
[16,] 0.5597850 0.88043000 0.44021500
[17,] 0.4976080 0.99521600 0.50239200
[18,] 0.4760361 0.95207212 0.52396394
[19,] 0.4884018 0.97680351 0.51159825
[20,] 0.4435946 0.88718924 0.55640538
[21,] 0.4149984 0.82999676 0.58500162
[22,] 0.3628367 0.72567346 0.63716327
[23,] 0.8211786 0.35764271 0.17882135
[24,] 0.8865671 0.22686571 0.11343286
[25,] 0.9022624 0.19547517 0.09773759
[26,] 0.9059394 0.18812118 0.09406059
[27,] 0.9537111 0.09257772 0.04628886
[28,] 0.9330101 0.13397975 0.06698987
[29,] 0.9095486 0.18090286 0.09045143
[30,] 0.8911606 0.21767873 0.10883937
[31,] 0.8557025 0.28859506 0.14429753
[32,] 0.8700308 0.25993839 0.12996919
[33,] 0.8187880 0.36242408 0.18121204
[34,] 0.9207933 0.15841347 0.07920673
[35,] 0.9757864 0.04842727 0.02421364
[36,] 0.9649294 0.07014117 0.03507058
[37,] 0.9517675 0.09646493 0.04823247
[38,] 0.9231081 0.15378377 0.07689189
[39,] 0.8790126 0.24197477 0.12098738
[40,] 0.8315531 0.33689372 0.16844686
[41,] 0.7611417 0.47771658 0.23885829
[42,] 0.8594436 0.28111284 0.14055642
[43,] 0.7725222 0.45495560 0.22747780
[44,] 0.8021928 0.39561433 0.19780717
[45,] 0.9256387 0.14872264 0.07436132
> postscript(file="/var/www/html/freestat/rcomp/tmp/104a91296726845.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/2qyqu1296726845.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/3ryjv1296726845.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/4xj7n1296726845.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/5exb91296726845.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 7
53.26393 169.80620 153.92940 123.64946 71.10994 123.58917 -165.73842
8 9 10 11 12 13 14
-147.63447 -19.04360 13.21707 -45.32772 -167.98657 -20.57352 54.59327
15 16 17 18 19 20 21
178.85068 258.54077 220.92954 248.68481 -74.83305 163.65688 61.47861
22 23 24 25 26 27 28
-91.28338 -45.94340 -84.91436 92.15158 129.08420 64.97943 179.10311
29 30 31 32 33 34 35
74.51704 331.83996 100.12407 -87.79510 -147.97670 12.20567 16.09569
36 37 38 39 40 41 42
-38.54232 -229.57621 -125.68877 -120.82346 93.92531 191.70028 -240.18008
43 44 45 46 47 48 49
-35.16066 -138.80898 -37.34822 -45.10816 25.29084 -118.36322 -161.69193
50 51 52 53 54 55 56
-27.68838 -59.87570 56.36503 49.41569 208.11816 -117.09199 -155.13363
57 58 59 60
-241.95199 -192.88039 -158.27433 -176.97704
> postscript(file="/var/www/html/freestat/rcomp/tmp/62xyj1296726845.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 53.26393 NA
1 169.80620 53.26393
2 153.92940 169.80620
3 123.64946 153.92940
4 71.10994 123.64946
5 123.58917 71.10994
6 -165.73842 123.58917
7 -147.63447 -165.73842
8 -19.04360 -147.63447
9 13.21707 -19.04360
10 -45.32772 13.21707
11 -167.98657 -45.32772
12 -20.57352 -167.98657
13 54.59327 -20.57352
14 178.85068 54.59327
15 258.54077 178.85068
16 220.92954 258.54077
17 248.68481 220.92954
18 -74.83305 248.68481
19 163.65688 -74.83305
20 61.47861 163.65688
21 -91.28338 61.47861
22 -45.94340 -91.28338
23 -84.91436 -45.94340
24 92.15158 -84.91436
25 129.08420 92.15158
26 64.97943 129.08420
27 179.10311 64.97943
28 74.51704 179.10311
29 331.83996 74.51704
30 100.12407 331.83996
31 -87.79510 100.12407
32 -147.97670 -87.79510
33 12.20567 -147.97670
34 16.09569 12.20567
35 -38.54232 16.09569
36 -229.57621 -38.54232
37 -125.68877 -229.57621
38 -120.82346 -125.68877
39 93.92531 -120.82346
40 191.70028 93.92531
41 -240.18008 191.70028
42 -35.16066 -240.18008
43 -138.80898 -35.16066
44 -37.34822 -138.80898
45 -45.10816 -37.34822
46 25.29084 -45.10816
47 -118.36322 25.29084
48 -161.69193 -118.36322
49 -27.68838 -161.69193
50 -59.87570 -27.68838
51 56.36503 -59.87570
52 49.41569 56.36503
53 208.11816 49.41569
54 -117.09199 208.11816
55 -155.13363 -117.09199
56 -241.95199 -155.13363
57 -192.88039 -241.95199
58 -158.27433 -192.88039
59 -176.97704 -158.27433
60 NA -176.97704
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 169.80620 53.26393
[2,] 153.92940 169.80620
[3,] 123.64946 153.92940
[4,] 71.10994 123.64946
[5,] 123.58917 71.10994
[6,] -165.73842 123.58917
[7,] -147.63447 -165.73842
[8,] -19.04360 -147.63447
[9,] 13.21707 -19.04360
[10,] -45.32772 13.21707
[11,] -167.98657 -45.32772
[12,] -20.57352 -167.98657
[13,] 54.59327 -20.57352
[14,] 178.85068 54.59327
[15,] 258.54077 178.85068
[16,] 220.92954 258.54077
[17,] 248.68481 220.92954
[18,] -74.83305 248.68481
[19,] 163.65688 -74.83305
[20,] 61.47861 163.65688
[21,] -91.28338 61.47861
[22,] -45.94340 -91.28338
[23,] -84.91436 -45.94340
[24,] 92.15158 -84.91436
[25,] 129.08420 92.15158
[26,] 64.97943 129.08420
[27,] 179.10311 64.97943
[28,] 74.51704 179.10311
[29,] 331.83996 74.51704
[30,] 100.12407 331.83996
[31,] -87.79510 100.12407
[32,] -147.97670 -87.79510
[33,] 12.20567 -147.97670
[34,] 16.09569 12.20567
[35,] -38.54232 16.09569
[36,] -229.57621 -38.54232
[37,] -125.68877 -229.57621
[38,] -120.82346 -125.68877
[39,] 93.92531 -120.82346
[40,] 191.70028 93.92531
[41,] -240.18008 191.70028
[42,] -35.16066 -240.18008
[43,] -138.80898 -35.16066
[44,] -37.34822 -138.80898
[45,] -45.10816 -37.34822
[46,] 25.29084 -45.10816
[47,] -118.36322 25.29084
[48,] -161.69193 -118.36322
[49,] -27.68838 -161.69193
[50,] -59.87570 -27.68838
[51,] 56.36503 -59.87570
[52,] 49.41569 56.36503
[53,] 208.11816 49.41569
[54,] -117.09199 208.11816
[55,] -155.13363 -117.09199
[56,] -241.95199 -155.13363
[57,] -192.88039 -241.95199
[58,] -158.27433 -192.88039
[59,] -176.97704 -158.27433
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 169.80620 53.26393
2 153.92940 169.80620
3 123.64946 153.92940
4 71.10994 123.64946
5 123.58917 71.10994
6 -165.73842 123.58917
7 -147.63447 -165.73842
8 -19.04360 -147.63447
9 13.21707 -19.04360
10 -45.32772 13.21707
11 -167.98657 -45.32772
12 -20.57352 -167.98657
13 54.59327 -20.57352
14 178.85068 54.59327
15 258.54077 178.85068
16 220.92954 258.54077
17 248.68481 220.92954
18 -74.83305 248.68481
19 163.65688 -74.83305
20 61.47861 163.65688
21 -91.28338 61.47861
22 -45.94340 -91.28338
23 -84.91436 -45.94340
24 92.15158 -84.91436
25 129.08420 92.15158
26 64.97943 129.08420
27 179.10311 64.97943
28 74.51704 179.10311
29 331.83996 74.51704
30 100.12407 331.83996
31 -87.79510 100.12407
32 -147.97670 -87.79510
33 12.20567 -147.97670
34 16.09569 12.20567
35 -38.54232 16.09569
36 -229.57621 -38.54232
37 -125.68877 -229.57621
38 -120.82346 -125.68877
39 93.92531 -120.82346
40 191.70028 93.92531
41 -240.18008 191.70028
42 -35.16066 -240.18008
43 -138.80898 -35.16066
44 -37.34822 -138.80898
45 -45.10816 -37.34822
46 25.29084 -45.10816
47 -118.36322 25.29084
48 -161.69193 -118.36322
49 -27.68838 -161.69193
50 -59.87570 -27.68838
51 56.36503 -59.87570
52 49.41569 56.36503
53 208.11816 49.41569
54 -117.09199 208.11816
55 -155.13363 -117.09199
56 -241.95199 -155.13363
57 -192.88039 -241.95199
58 -158.27433 -192.88039
59 -176.97704 -158.27433
> 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/7dsz01296726845.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/8uwwe1296726845.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/92s171296726845.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/10k5du1296726845.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/11cnd41296726845.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/12he1p1296726845.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/13kfgn1296726845.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/14ppzk1296726845.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/15ypqx1296726845.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/16h26r1296726845.tab")
+ }
>
> try(system("convert tmp/104a91296726845.ps tmp/104a91296726845.png",intern=TRUE))
character(0)
> try(system("convert tmp/2qyqu1296726845.ps tmp/2qyqu1296726845.png",intern=TRUE))
character(0)
> try(system("convert tmp/3ryjv1296726845.ps tmp/3ryjv1296726845.png",intern=TRUE))
character(0)
> try(system("convert tmp/4xj7n1296726845.ps tmp/4xj7n1296726845.png",intern=TRUE))
character(0)
> try(system("convert tmp/5exb91296726845.ps tmp/5exb91296726845.png",intern=TRUE))
character(0)
> try(system("convert tmp/62xyj1296726845.ps tmp/62xyj1296726845.png",intern=TRUE))
character(0)
> try(system("convert tmp/7dsz01296726845.ps tmp/7dsz01296726845.png",intern=TRUE))
character(0)
> try(system("convert tmp/8uwwe1296726845.ps tmp/8uwwe1296726845.png",intern=TRUE))
character(0)
> try(system("convert tmp/92s171296726845.ps tmp/92s171296726845.png",intern=TRUE))
character(0)
> try(system("convert tmp/10k5du1296726845.ps tmp/10k5du1296726845.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.946 2.531 4.496