R version 2.10.1 (2009-12-14)
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(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 = '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
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
1 0 0 0 0 0 0 0
2 0 0 0 0 0 0 0
3 0 0 0 0 0 0 0
4 0 0 0 0 0 0 0
5 1 0 0 0 0 0 0
6 0 1 0 0 0 0 0
7 0 0 1 0 0 0 0
8 0 0 0 1 0 0 0
9 0 0 0 0 1 0 0
10 0 0 0 0 0 1 0
11 0 0 0 0 0 0 1
12 0 0 0 0 0 0 0
13 0 0 0 0 0 0 0
14 0 0 0 0 0 0 0
15 0 0 0 0 0 0 0
16 0 0 0 0 0 0 0
17 1 0 0 0 0 0 0
18 0 1 0 0 0 0 0
19 0 0 1 0 0 0 0
20 0 0 0 1 0 0 0
21 0 0 0 0 1 0 0
22 0 0 0 0 0 1 0
23 0 0 0 0 0 0 1
24 0 0 0 0 0 0 0
25 0 0 0 0 0 0 0
26 0 0 0 0 0 0 0
27 0 0 0 0 0 0 0
28 0 0 0 0 0 0 0
29 1 0 0 0 0 0 0
30 0 1 0 0 0 0 0
31 0 0 1 0 0 0 0
32 0 0 0 1 0 0 0
33 0 0 0 0 1 0 0
34 0 0 0 0 0 1 0
35 0 0 0 0 0 0 1
36 0 0 0 0 0 0 0
37 0 0 0 0 0 0 0
38 0 0 0 0 0 0 0
39 0 0 0 0 0 0 0
40 0 0 0 0 0 0 0
41 1 0 0 0 0 0 0
42 0 1 0 0 0 0 0
43 0 0 1 0 0 0 0
44 0 0 0 1 0 0 0
45 0 0 0 0 1 0 0
46 0 0 0 0 0 1 0
47 0 0 0 0 0 0 1
48 0 0 0 0 0 0 0
49 0 0 0 0 0 0 0
50 0 0 0 0 0 0 0
51 0 0 0 0 0 0 0
52 0 0 0 0 0 0 0
53 1 0 0 0 0 0 0
54 0 1 0 0 0 0 0
55 0 0 1 0 0 0 0
56 0 0 0 1 0 0 0
57 0 0 0 0 1 0 0
58 0 0 0 0 0 1 0
59 0 0 0 0 0 0 1
60 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) Bewoonbare_opp Rentevoet Inflatie Werkloosheid
-3.687e+02 1.085e-02 -1.894e+01 -5.395e+00 -4.250e-04
M1 M2 M3 M4 M5
2.382e+01 1.168e+02 1.510e+02 3.384e+02 2.908e+02
M6 M7 M8 M9 M10
2.780e+02 3.653e+01 2.259e+01 3.260e-01 2.706e+01
M11
1.975e+02
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-356.834 -57.642 -6.291 69.286 189.633
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -3.687e+02 1.369e+03 -0.269 0.788987
Bewoonbare_opp 1.085e-02 3.449e-04 31.472 < 2e-16 ***
Rentevoet -1.894e+01 1.201e+02 -0.158 0.875444
Inflatie -5.395e+00 1.138e+01 -0.474 0.637871
Werkloosheid -4.250e-04 1.374e-03 -0.309 0.758553
M1 2.382e+01 7.300e+01 0.326 0.745742
M2 1.168e+02 7.793e+01 1.499 0.141021
M3 1.510e+02 8.845e+01 1.707 0.094856 .
M4 3.384e+02 9.557e+01 3.540 0.000958 ***
M5 2.908e+02 9.219e+01 3.155 0.002896 **
M6 2.780e+02 9.570e+01 2.905 0.005723 **
M7 3.653e+01 1.062e+02 0.344 0.732464
M8 2.259e+01 1.130e+02 0.200 0.842444
M9 3.260e-01 1.281e+02 0.003 0.997982
M10 2.706e+01 1.253e+02 0.216 0.829944
M11 1.975e+02 8.772e+01 2.251 0.029427 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 113.3 on 44 degrees of freedom
Multiple R-squared: 0.9812, Adjusted R-squared: 0.9747
F-statistic: 152.8 on 15 and 44 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.06374921 0.1274984 0.9362508
[2,] 0.05300232 0.1060046 0.9469977
[3,] 0.25183491 0.5036698 0.7481651
[4,] 0.46277660 0.9255532 0.5372234
[5,] 0.35858006 0.7171601 0.6414199
[6,] 0.25853452 0.5170690 0.7414655
[7,] 0.22443686 0.4488737 0.7755631
[8,] 0.17082273 0.3416455 0.8291773
[9,] 0.14299119 0.2859824 0.8570088
[10,] 0.09320205 0.1864041 0.9067979
[11,] 0.05913556 0.1182711 0.9408644
[12,] 0.21610999 0.4322200 0.7838900
[13,] 0.35097041 0.7019408 0.6490296
[14,] 0.32323045 0.6464609 0.6767695
[15,] 0.29963791 0.5992758 0.7003621
[16,] 0.45565426 0.9113085 0.5443457
[17,] 0.40609590 0.8121918 0.5939041
[18,] 0.30763910 0.6152782 0.6923609
[19,] 0.49430527 0.9886105 0.5056947
[20,] 0.46134107 0.9226821 0.5386589
[21,] 0.40293595 0.8058719 0.5970640
[22,] 0.30249249 0.6049850 0.6975075
[23,] 0.86193581 0.2761284 0.1380642
> postscript(file="/var/yougetitorg/rcomp/tmp/15l4v1296730703.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/yougetitorg/rcomp/tmp/229bn1296730703.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/yougetitorg/rcomp/tmp/3nssv1296730703.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/yougetitorg/rcomp/tmp/4hpnw1296730703.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/yougetitorg/rcomp/tmp/5r6am1296730703.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
98.135386 137.387860 110.286524 -13.055650 -76.847082 -27.682933
7 8 9 10 11 12
-202.503553 -92.287137 17.688622 45.380397 -6.247955 -95.489748
13 14 15 16 17 18
11.615693 -39.974451 75.045280 92.807239 49.793668 73.191500
19 20 21 22 23 24
-52.299794 189.633320 49.442387 -92.368497 -64.283894 -21.122178
25 26 27 28 29 30
113.611505 52.437584 -2.844370 7.456172 -30.712391 143.314512
31 32 33 34 35 36
156.175292 -35.039270 -68.116041 67.983689 91.948439 148.254776
37 38 39 40 41 42
-160.711936 -103.887762 -135.983877 -33.716385 84.893667 -356.834357
43 44 45 46 47 48
51.589167 -55.972006 97.906780 58.587896 49.082610 -3.452863
49 50 51 52 53 54
-62.650647 -45.963231 -46.503558 -53.491376 -27.127861 168.011278
55 56 57 58 59 60
47.038888 -6.334907 -96.921748 -79.583485 -70.499199 -28.189987
> postscript(file="/var/yougetitorg/rcomp/tmp/62und1296730703.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 98.135386 NA
1 137.387860 98.135386
2 110.286524 137.387860
3 -13.055650 110.286524
4 -76.847082 -13.055650
5 -27.682933 -76.847082
6 -202.503553 -27.682933
7 -92.287137 -202.503553
8 17.688622 -92.287137
9 45.380397 17.688622
10 -6.247955 45.380397
11 -95.489748 -6.247955
12 11.615693 -95.489748
13 -39.974451 11.615693
14 75.045280 -39.974451
15 92.807239 75.045280
16 49.793668 92.807239
17 73.191500 49.793668
18 -52.299794 73.191500
19 189.633320 -52.299794
20 49.442387 189.633320
21 -92.368497 49.442387
22 -64.283894 -92.368497
23 -21.122178 -64.283894
24 113.611505 -21.122178
25 52.437584 113.611505
26 -2.844370 52.437584
27 7.456172 -2.844370
28 -30.712391 7.456172
29 143.314512 -30.712391
30 156.175292 143.314512
31 -35.039270 156.175292
32 -68.116041 -35.039270
33 67.983689 -68.116041
34 91.948439 67.983689
35 148.254776 91.948439
36 -160.711936 148.254776
37 -103.887762 -160.711936
38 -135.983877 -103.887762
39 -33.716385 -135.983877
40 84.893667 -33.716385
41 -356.834357 84.893667
42 51.589167 -356.834357
43 -55.972006 51.589167
44 97.906780 -55.972006
45 58.587896 97.906780
46 49.082610 58.587896
47 -3.452863 49.082610
48 -62.650647 -3.452863
49 -45.963231 -62.650647
50 -46.503558 -45.963231
51 -53.491376 -46.503558
52 -27.127861 -53.491376
53 168.011278 -27.127861
54 47.038888 168.011278
55 -6.334907 47.038888
56 -96.921748 -6.334907
57 -79.583485 -96.921748
58 -70.499199 -79.583485
59 -28.189987 -70.499199
60 NA -28.189987
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 137.387860 98.135386
[2,] 110.286524 137.387860
[3,] -13.055650 110.286524
[4,] -76.847082 -13.055650
[5,] -27.682933 -76.847082
[6,] -202.503553 -27.682933
[7,] -92.287137 -202.503553
[8,] 17.688622 -92.287137
[9,] 45.380397 17.688622
[10,] -6.247955 45.380397
[11,] -95.489748 -6.247955
[12,] 11.615693 -95.489748
[13,] -39.974451 11.615693
[14,] 75.045280 -39.974451
[15,] 92.807239 75.045280
[16,] 49.793668 92.807239
[17,] 73.191500 49.793668
[18,] -52.299794 73.191500
[19,] 189.633320 -52.299794
[20,] 49.442387 189.633320
[21,] -92.368497 49.442387
[22,] -64.283894 -92.368497
[23,] -21.122178 -64.283894
[24,] 113.611505 -21.122178
[25,] 52.437584 113.611505
[26,] -2.844370 52.437584
[27,] 7.456172 -2.844370
[28,] -30.712391 7.456172
[29,] 143.314512 -30.712391
[30,] 156.175292 143.314512
[31,] -35.039270 156.175292
[32,] -68.116041 -35.039270
[33,] 67.983689 -68.116041
[34,] 91.948439 67.983689
[35,] 148.254776 91.948439
[36,] -160.711936 148.254776
[37,] -103.887762 -160.711936
[38,] -135.983877 -103.887762
[39,] -33.716385 -135.983877
[40,] 84.893667 -33.716385
[41,] -356.834357 84.893667
[42,] 51.589167 -356.834357
[43,] -55.972006 51.589167
[44,] 97.906780 -55.972006
[45,] 58.587896 97.906780
[46,] 49.082610 58.587896
[47,] -3.452863 49.082610
[48,] -62.650647 -3.452863
[49,] -45.963231 -62.650647
[50,] -46.503558 -45.963231
[51,] -53.491376 -46.503558
[52,] -27.127861 -53.491376
[53,] 168.011278 -27.127861
[54,] 47.038888 168.011278
[55,] -6.334907 47.038888
[56,] -96.921748 -6.334907
[57,] -79.583485 -96.921748
[58,] -70.499199 -79.583485
[59,] -28.189987 -70.499199
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 137.387860 98.135386
2 110.286524 137.387860
3 -13.055650 110.286524
4 -76.847082 -13.055650
5 -27.682933 -76.847082
6 -202.503553 -27.682933
7 -92.287137 -202.503553
8 17.688622 -92.287137
9 45.380397 17.688622
10 -6.247955 45.380397
11 -95.489748 -6.247955
12 11.615693 -95.489748
13 -39.974451 11.615693
14 75.045280 -39.974451
15 92.807239 75.045280
16 49.793668 92.807239
17 73.191500 49.793668
18 -52.299794 73.191500
19 189.633320 -52.299794
20 49.442387 189.633320
21 -92.368497 49.442387
22 -64.283894 -92.368497
23 -21.122178 -64.283894
24 113.611505 -21.122178
25 52.437584 113.611505
26 -2.844370 52.437584
27 7.456172 -2.844370
28 -30.712391 7.456172
29 143.314512 -30.712391
30 156.175292 143.314512
31 -35.039270 156.175292
32 -68.116041 -35.039270
33 67.983689 -68.116041
34 91.948439 67.983689
35 148.254776 91.948439
36 -160.711936 148.254776
37 -103.887762 -160.711936
38 -135.983877 -103.887762
39 -33.716385 -135.983877
40 84.893667 -33.716385
41 -356.834357 84.893667
42 51.589167 -356.834357
43 -55.972006 51.589167
44 97.906780 -55.972006
45 58.587896 97.906780
46 49.082610 58.587896
47 -3.452863 49.082610
48 -62.650647 -3.452863
49 -45.963231 -62.650647
50 -46.503558 -45.963231
51 -53.491376 -46.503558
52 -27.127861 -53.491376
53 168.011278 -27.127861
54 47.038888 168.011278
55 -6.334907 47.038888
56 -96.921748 -6.334907
57 -79.583485 -96.921748
58 -70.499199 -79.583485
59 -28.189987 -70.499199
> 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/yougetitorg/rcomp/tmp/7cat41296730703.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/yougetitorg/rcomp/tmp/8dx901296730703.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/yougetitorg/rcomp/tmp/958u81296730703.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/yougetitorg/rcomp/tmp/10xoj01296730703.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/yougetitorg/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/yougetitorg/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/yougetitorg/rcomp/tmp/11ptpm1296730703.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/yougetitorg/rcomp/tmp/121kyv1296730703.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/yougetitorg/rcomp/tmp/13naia1296730703.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/yougetitorg/rcomp/tmp/14og2z1296730703.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/yougetitorg/rcomp/tmp/15vpxl1296730703.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/yougetitorg/rcomp/tmp/166m3h1296730703.tab")
+ }
>
> try(system("convert tmp/15l4v1296730703.ps tmp/15l4v1296730703.png",intern=TRUE))
character(0)
> try(system("convert tmp/229bn1296730703.ps tmp/229bn1296730703.png",intern=TRUE))
character(0)
> try(system("convert tmp/3nssv1296730703.ps tmp/3nssv1296730703.png",intern=TRUE))
character(0)
> try(system("convert tmp/4hpnw1296730703.ps tmp/4hpnw1296730703.png",intern=TRUE))
character(0)
> try(system("convert tmp/5r6am1296730703.ps tmp/5r6am1296730703.png",intern=TRUE))
character(0)
> try(system("convert tmp/62und1296730703.ps tmp/62und1296730703.png",intern=TRUE))
character(0)
> try(system("convert tmp/7cat41296730703.ps tmp/7cat41296730703.png",intern=TRUE))
character(0)
> try(system("convert tmp/8dx901296730703.ps tmp/8dx901296730703.png",intern=TRUE))
character(0)
> try(system("convert tmp/958u81296730703.ps tmp/958u81296730703.png",intern=TRUE))
character(0)
> try(system("convert tmp/10xoj01296730703.ps tmp/10xoj01296730703.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.300 2.130 4.447