R version 2.9.0 (2009-04-17)
Copyright (C) 2009 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> x <- array(list(101.6
+ ,595454
+ ,594167
+ ,611324
+ ,612613
+ ,610763
+ ,94.6
+ ,590865
+ ,595454
+ ,594167
+ ,611324
+ ,612613
+ ,95.9
+ ,589379
+ ,590865
+ ,595454
+ ,594167
+ ,611324
+ ,104.7
+ ,584428
+ ,589379
+ ,590865
+ ,595454
+ ,594167
+ ,102.8
+ ,573100
+ ,584428
+ ,589379
+ ,590865
+ ,595454
+ ,98.1
+ ,567456
+ ,573100
+ ,584428
+ ,589379
+ ,590865
+ ,113.9
+ ,569028
+ ,567456
+ ,573100
+ ,584428
+ ,589379
+ ,80.9
+ ,620735
+ ,569028
+ ,567456
+ ,573100
+ ,584428
+ ,95.7
+ ,628884
+ ,620735
+ ,569028
+ ,567456
+ ,573100
+ ,113.2
+ ,628232
+ ,628884
+ ,620735
+ ,569028
+ ,567456
+ ,105.9
+ ,612117
+ ,628232
+ ,628884
+ ,620735
+ ,569028
+ ,108.8
+ ,595404
+ ,612117
+ ,628232
+ ,628884
+ ,620735
+ ,102.3
+ ,597141
+ ,595404
+ ,612117
+ ,628232
+ ,628884
+ ,99
+ ,593408
+ ,597141
+ ,595404
+ ,612117
+ ,628232
+ ,100.7
+ ,590072
+ ,593408
+ ,597141
+ ,595404
+ ,612117
+ ,115.5
+ ,579799
+ ,590072
+ ,593408
+ ,597141
+ ,595404
+ ,100.7
+ ,574205
+ ,579799
+ ,590072
+ ,593408
+ ,597141
+ ,109.9
+ ,572775
+ ,574205
+ ,579799
+ ,590072
+ ,593408
+ ,114.6
+ ,572942
+ ,572775
+ ,574205
+ ,579799
+ ,590072
+ ,85.4
+ ,619567
+ ,572942
+ ,572775
+ ,574205
+ ,579799
+ ,100.5
+ ,625809
+ ,619567
+ ,572942
+ ,572775
+ ,574205
+ ,114.8
+ ,619916
+ ,625809
+ ,619567
+ ,572942
+ ,572775
+ ,116.5
+ ,587625
+ ,619916
+ ,625809
+ ,619567
+ ,572942
+ ,112.9
+ ,565742
+ ,587625
+ ,619916
+ ,625809
+ ,619567
+ ,102
+ ,557274
+ ,565742
+ ,587625
+ ,619916
+ ,625809
+ ,106
+ ,560576
+ ,557274
+ ,565742
+ ,587625
+ ,619916
+ ,105.3
+ ,548854
+ ,560576
+ ,557274
+ ,565742
+ ,587625
+ ,118.8
+ ,531673
+ ,548854
+ ,560576
+ ,557274
+ ,565742
+ ,106.1
+ ,525919
+ ,531673
+ ,548854
+ ,560576
+ ,557274
+ ,109.3
+ ,511038
+ ,525919
+ ,531673
+ ,548854
+ ,560576
+ ,117.2
+ ,498662
+ ,511038
+ ,525919
+ ,531673
+ ,548854
+ ,92.5
+ ,555362
+ ,498662
+ ,511038
+ ,525919
+ ,531673
+ ,104.2
+ ,564591
+ ,555362
+ ,498662
+ ,511038
+ ,525919
+ ,112.5
+ ,541657
+ ,564591
+ ,555362
+ ,498662
+ ,511038
+ ,122.4
+ ,527070
+ ,541657
+ ,564591
+ ,555362
+ ,498662
+ ,113.3
+ ,509846
+ ,527070
+ ,541657
+ ,564591
+ ,555362
+ ,100
+ ,514258
+ ,509846
+ ,527070
+ ,541657
+ ,564591
+ ,110.7
+ ,516922
+ ,514258
+ ,509846
+ ,527070
+ ,541657
+ ,112.8
+ ,507561
+ ,516922
+ ,514258
+ ,509846
+ ,527070
+ ,109.8
+ ,492622
+ ,507561
+ ,516922
+ ,514258
+ ,509846
+ ,117.3
+ ,490243
+ ,492622
+ ,507561
+ ,516922
+ ,514258
+ ,109.1
+ ,469357
+ ,490243
+ ,492622
+ ,507561
+ ,516922
+ ,115.9
+ ,477580
+ ,469357
+ ,490243
+ ,492622
+ ,507561
+ ,96
+ ,528379
+ ,477580
+ ,469357
+ ,490243
+ ,492622
+ ,99.8
+ ,533590
+ ,528379
+ ,477580
+ ,469357
+ ,490243
+ ,116.8
+ ,517945
+ ,533590
+ ,528379
+ ,477580
+ ,469357
+ ,115.7
+ ,506174
+ ,517945
+ ,533590
+ ,528379
+ ,477580
+ ,99.4
+ ,501866
+ ,506174
+ ,517945
+ ,533590
+ ,528379
+ ,94.3
+ ,516141
+ ,501866
+ ,506174
+ ,517945
+ ,533590
+ ,91
+ ,528222
+ ,516141
+ ,501866
+ ,506174
+ ,517945
+ ,93.2
+ ,532638
+ ,528222
+ ,516141
+ ,501866
+ ,506174
+ ,103.1
+ ,536322
+ ,532638
+ ,528222
+ ,516141
+ ,501866
+ ,94.1
+ ,536535
+ ,536322
+ ,532638
+ ,528222
+ ,516141
+ ,91.8
+ ,523597
+ ,536535
+ ,536322
+ ,532638
+ ,528222
+ ,102.7
+ ,536214
+ ,523597
+ ,536535
+ ,536322
+ ,532638
+ ,82.6
+ ,586570
+ ,536214
+ ,523597
+ ,536535
+ ,536322
+ ,89.1
+ ,596594
+ ,586570
+ ,536214
+ ,523597
+ ,536535)
+ ,dim=c(6
+ ,57)
+ ,dimnames=list(c('X'
+ ,'Niet_werkende_werkzoekenden'
+ ,'y1'
+ ,'y2'
+ ,'y3'
+ ,'y4')
+ ,1:57))
> y <- array(NA,dim=c(6,57),dimnames=list(c('X','Niet_werkende_werkzoekenden','y1','y2','y3','y4'),1:57))
> 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 = '2'
> #'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
Niet_werkende_werkzoekenden X y1 y2 y3 y4 M1 M2 M3 M4 M5
1 595454 101.6 594167 611324 612613 610763 1 0 0 0 0
2 590865 94.6 595454 594167 611324 612613 0 1 0 0 0
3 589379 95.9 590865 595454 594167 611324 0 0 1 0 0
4 584428 104.7 589379 590865 595454 594167 0 0 0 1 0
5 573100 102.8 584428 589379 590865 595454 0 0 0 0 1
6 567456 98.1 573100 584428 589379 590865 0 0 0 0 0
7 569028 113.9 567456 573100 584428 589379 0 0 0 0 0
8 620735 80.9 569028 567456 573100 584428 0 0 0 0 0
9 628884 95.7 620735 569028 567456 573100 0 0 0 0 0
10 628232 113.2 628884 620735 569028 567456 0 0 0 0 0
11 612117 105.9 628232 628884 620735 569028 0 0 0 0 0
12 595404 108.8 612117 628232 628884 620735 0 0 0 0 0
13 597141 102.3 595404 612117 628232 628884 1 0 0 0 0
14 593408 99.0 597141 595404 612117 628232 0 1 0 0 0
15 590072 100.7 593408 597141 595404 612117 0 0 1 0 0
16 579799 115.5 590072 593408 597141 595404 0 0 0 1 0
17 574205 100.7 579799 590072 593408 597141 0 0 0 0 1
18 572775 109.9 574205 579799 590072 593408 0 0 0 0 0
19 572942 114.6 572775 574205 579799 590072 0 0 0 0 0
20 619567 85.4 572942 572775 574205 579799 0 0 0 0 0
21 625809 100.5 619567 572942 572775 574205 0 0 0 0 0
22 619916 114.8 625809 619567 572942 572775 0 0 0 0 0
23 587625 116.5 619916 625809 619567 572942 0 0 0 0 0
24 565742 112.9 587625 619916 625809 619567 0 0 0 0 0
25 557274 102.0 565742 587625 619916 625809 1 0 0 0 0
26 560576 106.0 557274 565742 587625 619916 0 1 0 0 0
27 548854 105.3 560576 557274 565742 587625 0 0 1 0 0
28 531673 118.8 548854 560576 557274 565742 0 0 0 1 0
29 525919 106.1 531673 548854 560576 557274 0 0 0 0 1
30 511038 109.3 525919 531673 548854 560576 0 0 0 0 0
31 498662 117.2 511038 525919 531673 548854 0 0 0 0 0
32 555362 92.5 498662 511038 525919 531673 0 0 0 0 0
33 564591 104.2 555362 498662 511038 525919 0 0 0 0 0
34 541657 112.5 564591 555362 498662 511038 0 0 0 0 0
35 527070 122.4 541657 564591 555362 498662 0 0 0 0 0
36 509846 113.3 527070 541657 564591 555362 0 0 0 0 0
37 514258 100.0 509846 527070 541657 564591 1 0 0 0 0
38 516922 110.7 514258 509846 527070 541657 0 1 0 0 0
39 507561 112.8 516922 514258 509846 527070 0 0 1 0 0
40 492622 109.8 507561 516922 514258 509846 0 0 0 1 0
41 490243 117.3 492622 507561 516922 514258 0 0 0 0 1
42 469357 109.1 490243 492622 507561 516922 0 0 0 0 0
43 477580 115.9 469357 490243 492622 507561 0 0 0 0 0
44 528379 96.0 477580 469357 490243 492622 0 0 0 0 0
45 533590 99.8 528379 477580 469357 490243 0 0 0 0 0
46 517945 116.8 533590 528379 477580 469357 0 0 0 0 0
47 506174 115.7 517945 533590 528379 477580 0 0 0 0 0
48 501866 99.4 506174 517945 533590 528379 0 0 0 0 0
49 516141 94.3 501866 506174 517945 533590 1 0 0 0 0
50 528222 91.0 516141 501866 506174 517945 0 1 0 0 0
51 532638 93.2 528222 516141 501866 506174 0 0 1 0 0
52 536322 103.1 532638 528222 516141 501866 0 0 0 1 0
53 536535 94.1 536322 532638 528222 516141 0 0 0 0 1
54 523597 91.8 536535 536322 532638 528222 0 0 0 0 0
55 536214 102.7 523597 536535 536322 532638 0 0 0 0 0
56 586570 82.6 536214 523597 536535 536322 0 0 0 0 0
57 596594 89.1 586570 536214 523597 536535 0 0 0 0 0
M6 M7 M8 M9 M10 M11 t
1 0 0 0 0 0 0 1
2 0 0 0 0 0 0 2
3 0 0 0 0 0 0 3
4 0 0 0 0 0 0 4
5 0 0 0 0 0 0 5
6 1 0 0 0 0 0 6
7 0 1 0 0 0 0 7
8 0 0 1 0 0 0 8
9 0 0 0 1 0 0 9
10 0 0 0 0 1 0 10
11 0 0 0 0 0 1 11
12 0 0 0 0 0 0 12
13 0 0 0 0 0 0 13
14 0 0 0 0 0 0 14
15 0 0 0 0 0 0 15
16 0 0 0 0 0 0 16
17 0 0 0 0 0 0 17
18 1 0 0 0 0 0 18
19 0 1 0 0 0 0 19
20 0 0 1 0 0 0 20
21 0 0 0 1 0 0 21
22 0 0 0 0 1 0 22
23 0 0 0 0 0 1 23
24 0 0 0 0 0 0 24
25 0 0 0 0 0 0 25
26 0 0 0 0 0 0 26
27 0 0 0 0 0 0 27
28 0 0 0 0 0 0 28
29 0 0 0 0 0 0 29
30 1 0 0 0 0 0 30
31 0 1 0 0 0 0 31
32 0 0 1 0 0 0 32
33 0 0 0 1 0 0 33
34 0 0 0 0 1 0 34
35 0 0 0 0 0 1 35
36 0 0 0 0 0 0 36
37 0 0 0 0 0 0 37
38 0 0 0 0 0 0 38
39 0 0 0 0 0 0 39
40 0 0 0 0 0 0 40
41 0 0 0 0 0 0 41
42 1 0 0 0 0 0 42
43 0 1 0 0 0 0 43
44 0 0 1 0 0 0 44
45 0 0 0 1 0 0 45
46 0 0 0 0 1 0 46
47 0 0 0 0 0 1 47
48 0 0 0 0 0 0 48
49 0 0 0 0 0 0 49
50 0 0 0 0 0 0 50
51 0 0 0 0 0 0 51
52 0 0 0 0 0 0 52
53 0 0 0 0 0 0 53
54 1 0 0 0 0 0 54
55 0 1 0 0 0 0 55
56 0 0 1 0 0 0 56
57 0 0 0 1 0 0 57
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X y1 y2 y3 y4
6.897e+04 -4.183e+02 9.305e-01 3.301e-02 1.879e-01 -2.198e-01
M1 M2 M3 M4 M5 M6
1.711e+04 1.821e+04 1.218e+04 7.222e+03 8.205e+03 3.015e+03
M7 M8 M9 M10 M11 t
2.021e+04 5.837e+04 2.379e+04 7.833e+03 -1.001e+04 -5.382e+01
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-11942.0 -2866.9 530.9 3563.2 13876.5
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.897e+04 4.253e+04 1.622 0.11293
X -4.183e+02 2.059e+02 -2.032 0.04905 *
y1 9.304e-01 1.622e-01 5.738 1.20e-06 ***
y2 3.301e-02 2.111e-01 0.156 0.87652
y3 1.879e-01 2.204e-01 0.853 0.39911
y4 -2.198e-01 1.437e-01 -1.529 0.13433
M1 1.711e+04 5.565e+03 3.075 0.00383 **
M2 1.821e+04 6.728e+03 2.706 0.01005 *
M3 1.218e+04 7.006e+03 1.738 0.09003 .
M4 7.222e+03 5.854e+03 1.234 0.22475
M5 8.205e+03 5.106e+03 1.607 0.11613
M6 3.015e+03 5.470e+03 0.551 0.58457
M7 2.021e+04 5.609e+03 3.603 0.00088 ***
M8 5.837e+04 7.108e+03 8.211 4.96e-10 ***
M9 2.379e+04 1.227e+04 1.939 0.05979 .
M10 7.833e+03 1.286e+04 0.609 0.54586
M11 -1.001e+04 8.790e+03 -1.139 0.26165
t -5.382e+01 1.129e+02 -0.477 0.63633
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 6734 on 39 degrees of freedom
Multiple R-squared: 0.982, Adjusted R-squared: 0.9741
F-statistic: 124.9 on 17 and 39 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.004328684 0.008657368 0.99567132
[2,] 0.136928882 0.273857765 0.86307112
[3,] 0.270614969 0.541229939 0.72938503
[4,] 0.187114406 0.374228812 0.81288559
[5,] 0.243700794 0.487401589 0.75629921
[6,] 0.388966722 0.777933444 0.61103328
[7,] 0.302034587 0.604069175 0.69796541
[8,] 0.239272177 0.478544355 0.76072782
[9,] 0.343585823 0.687171647 0.65641418
[10,] 0.420327385 0.840654771 0.57967261
[11,] 0.546924484 0.906151032 0.45307552
[12,] 0.764580774 0.470838452 0.23541923
[13,] 0.885226968 0.229546064 0.11477303
[14,] 0.914183531 0.171632937 0.08581647
[15,] 0.879853314 0.240293372 0.12014669
[16,] 0.767751501 0.464496997 0.23224850
> postscript(file="/var/www/html/rcomp/tmp/1pwh31258731432.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> plot(x[,1], type='l', main='Actuals and Interpolation', ylab='value of Actuals and Interpolation (dots)', xlab='time or index')
> points(x[,1]-mysum$resid)
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/23r791258731432.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> plot(mysum$resid, type='b', pch=19, main='Residuals', ylab='value of Residuals', xlab='time or index')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/3tj1g1258731432.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> hist(mysum$resid, main='Residual Histogram', xlab='values of Residuals')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/4vopj1258731432.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> densityplot(~mysum$resid,col='black',main='Residual Density Plot', xlab='values of Residuals')
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/5fimt1258731432.ps",horizontal=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 = 57
Frequency = 1
1 2 3 4 5 6
-1965.8760 -10503.7182 1803.4188 3066.0608 -4185.7580 3422.4248
7 8 9 10 11 12
695.2897 257.4087 -370.5241 11487.7391 1186.0303 577.6398
13 14 15 16 17 18
530.8613 -3801.0388 2670.1200 2826.6359 863.8585 13876.5096
19 20 21 22 23 24
1584.3996 -3424.8126 -590.7909 7820.2323 -9307.3719 -3339.3118
25 26 27 28 29 30
-9521.5711 7788.2752 -3923.2248 -2866.9299 -971.5444 -420.9568
31 32 33 34 35 36
-11941.9776 5632.6262 3563.2473 -11288.1070 3825.8382 -2104.7232
37 38 39 40 41 42
2527.8606 2792.6073 -2202.4898 -9377.6872 5128.9133 -8892.1408
43 44 45 46 47 48
5299.3510 -128.7248 -2840.0447 -8019.8644 4295.5034 4866.3952
49 50 51 52 53 54
8428.7252 3723.8745 1652.1758 6351.9203 -835.4694 -7985.8369
55 56 57
4362.9373 -2336.4975 238.1125
> postscript(file="/var/www/html/rcomp/tmp/69yzh1258731432.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> dum <- cbind(lag(myerror,k=1),myerror)
> dum
Time Series:
Start = 0
End = 57
Frequency = 1
lag(myerror, k = 1) myerror
0 -1965.8760 NA
1 -10503.7182 -1965.8760
2 1803.4188 -10503.7182
3 3066.0608 1803.4188
4 -4185.7580 3066.0608
5 3422.4248 -4185.7580
6 695.2897 3422.4248
7 257.4087 695.2897
8 -370.5241 257.4087
9 11487.7391 -370.5241
10 1186.0303 11487.7391
11 577.6398 1186.0303
12 530.8613 577.6398
13 -3801.0388 530.8613
14 2670.1200 -3801.0388
15 2826.6359 2670.1200
16 863.8585 2826.6359
17 13876.5096 863.8585
18 1584.3996 13876.5096
19 -3424.8126 1584.3996
20 -590.7909 -3424.8126
21 7820.2323 -590.7909
22 -9307.3719 7820.2323
23 -3339.3118 -9307.3719
24 -9521.5711 -3339.3118
25 7788.2752 -9521.5711
26 -3923.2248 7788.2752
27 -2866.9299 -3923.2248
28 -971.5444 -2866.9299
29 -420.9568 -971.5444
30 -11941.9776 -420.9568
31 5632.6262 -11941.9776
32 3563.2473 5632.6262
33 -11288.1070 3563.2473
34 3825.8382 -11288.1070
35 -2104.7232 3825.8382
36 2527.8606 -2104.7232
37 2792.6073 2527.8606
38 -2202.4898 2792.6073
39 -9377.6872 -2202.4898
40 5128.9133 -9377.6872
41 -8892.1408 5128.9133
42 5299.3510 -8892.1408
43 -128.7248 5299.3510
44 -2840.0447 -128.7248
45 -8019.8644 -2840.0447
46 4295.5034 -8019.8644
47 4866.3952 4295.5034
48 8428.7252 4866.3952
49 3723.8745 8428.7252
50 1652.1758 3723.8745
51 6351.9203 1652.1758
52 -835.4694 6351.9203
53 -7985.8369 -835.4694
54 4362.9373 -7985.8369
55 -2336.4975 4362.9373
56 238.1125 -2336.4975
57 NA 238.1125
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -10503.7182 -1965.8760
[2,] 1803.4188 -10503.7182
[3,] 3066.0608 1803.4188
[4,] -4185.7580 3066.0608
[5,] 3422.4248 -4185.7580
[6,] 695.2897 3422.4248
[7,] 257.4087 695.2897
[8,] -370.5241 257.4087
[9,] 11487.7391 -370.5241
[10,] 1186.0303 11487.7391
[11,] 577.6398 1186.0303
[12,] 530.8613 577.6398
[13,] -3801.0388 530.8613
[14,] 2670.1200 -3801.0388
[15,] 2826.6359 2670.1200
[16,] 863.8585 2826.6359
[17,] 13876.5096 863.8585
[18,] 1584.3996 13876.5096
[19,] -3424.8126 1584.3996
[20,] -590.7909 -3424.8126
[21,] 7820.2323 -590.7909
[22,] -9307.3719 7820.2323
[23,] -3339.3118 -9307.3719
[24,] -9521.5711 -3339.3118
[25,] 7788.2752 -9521.5711
[26,] -3923.2248 7788.2752
[27,] -2866.9299 -3923.2248
[28,] -971.5444 -2866.9299
[29,] -420.9568 -971.5444
[30,] -11941.9776 -420.9568
[31,] 5632.6262 -11941.9776
[32,] 3563.2473 5632.6262
[33,] -11288.1070 3563.2473
[34,] 3825.8382 -11288.1070
[35,] -2104.7232 3825.8382
[36,] 2527.8606 -2104.7232
[37,] 2792.6073 2527.8606
[38,] -2202.4898 2792.6073
[39,] -9377.6872 -2202.4898
[40,] 5128.9133 -9377.6872
[41,] -8892.1408 5128.9133
[42,] 5299.3510 -8892.1408
[43,] -128.7248 5299.3510
[44,] -2840.0447 -128.7248
[45,] -8019.8644 -2840.0447
[46,] 4295.5034 -8019.8644
[47,] 4866.3952 4295.5034
[48,] 8428.7252 4866.3952
[49,] 3723.8745 8428.7252
[50,] 1652.1758 3723.8745
[51,] 6351.9203 1652.1758
[52,] -835.4694 6351.9203
[53,] -7985.8369 -835.4694
[54,] 4362.9373 -7985.8369
[55,] -2336.4975 4362.9373
[56,] 238.1125 -2336.4975
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -10503.7182 -1965.8760
2 1803.4188 -10503.7182
3 3066.0608 1803.4188
4 -4185.7580 3066.0608
5 3422.4248 -4185.7580
6 695.2897 3422.4248
7 257.4087 695.2897
8 -370.5241 257.4087
9 11487.7391 -370.5241
10 1186.0303 11487.7391
11 577.6398 1186.0303
12 530.8613 577.6398
13 -3801.0388 530.8613
14 2670.1200 -3801.0388
15 2826.6359 2670.1200
16 863.8585 2826.6359
17 13876.5096 863.8585
18 1584.3996 13876.5096
19 -3424.8126 1584.3996
20 -590.7909 -3424.8126
21 7820.2323 -590.7909
22 -9307.3719 7820.2323
23 -3339.3118 -9307.3719
24 -9521.5711 -3339.3118
25 7788.2752 -9521.5711
26 -3923.2248 7788.2752
27 -2866.9299 -3923.2248
28 -971.5444 -2866.9299
29 -420.9568 -971.5444
30 -11941.9776 -420.9568
31 5632.6262 -11941.9776
32 3563.2473 5632.6262
33 -11288.1070 3563.2473
34 3825.8382 -11288.1070
35 -2104.7232 3825.8382
36 2527.8606 -2104.7232
37 2792.6073 2527.8606
38 -2202.4898 2792.6073
39 -9377.6872 -2202.4898
40 5128.9133 -9377.6872
41 -8892.1408 5128.9133
42 5299.3510 -8892.1408
43 -128.7248 5299.3510
44 -2840.0447 -128.7248
45 -8019.8644 -2840.0447
46 4295.5034 -8019.8644
47 4866.3952 4295.5034
48 8428.7252 4866.3952
49 3723.8745 8428.7252
50 1652.1758 3723.8745
51 6351.9203 1652.1758
52 -835.4694 6351.9203
53 -7985.8369 -835.4694
54 4362.9373 -7985.8369
55 -2336.4975 4362.9373
56 238.1125 -2336.4975
> plot(z,main=paste('Residual Lag plot, lowess, and regression line'), ylab='values of Residuals', xlab='lagged values of Residuals')
> lines(lowess(z))
> abline(lm(z))
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/7kbxq1258731432.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> acf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Autocorrelation Function')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/89zjc1258731432.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> pacf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Partial Autocorrelation Function')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/99zx01258731432.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> opar <- par(mfrow = c(2,2), oma = c(0, 0, 1.1, 0))
> plot(mylm, las = 1, sub='Residual Diagnostics')
> par(opar)
> dev.off()
null device
1
> if (n > n25) {
+ postscript(file="/var/www/html/rcomp/tmp/10hu6a1258731432.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
+ plot(kp3:nmkm3,gqarr[,2], main='Goldfeld-Quandt test',ylab='2-sided p-value',xlab='breakpoint')
+ grid()
+ dev.off()
+ }
null device
1
>
> #Note: the /var/www/html/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/html/rcomp/createtable")
>
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Estimated Regression Equation', 1, TRUE)
> a<-table.row.end(a)
> myeq <- colnames(x)[1]
> myeq <- paste(myeq, '[t] = ', sep='')
> for (i in 1:k){
+ if (mysum$coefficients[i,1] > 0) myeq <- paste(myeq, '+', '')
+ myeq <- paste(myeq, mysum$coefficients[i,1], sep=' ')
+ if (rownames(mysum$coefficients)[i] != '(Intercept)') {
+ myeq <- paste(myeq, rownames(mysum$coefficients)[i], sep='')
+ if (rownames(mysum$coefficients)[i] != 't') myeq <- paste(myeq, '[t]', sep='')
+ }
+ }
> myeq <- paste(myeq, ' + e[t]')
> a<-table.row.start(a)
> a<-table.element(a, myeq)
> a<-table.row.end(a)
> a<-table.end(a)
> table.save(a,file="/var/www/html/rcomp/tmp/1120kp1258731433.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a,hyperlink('http://www.xycoon.com/ols1.htm','Multiple Linear Regression - Ordinary Least Squares',''), 6, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a,'Variable',header=TRUE)
> a<-table.element(a,'Parameter',header=TRUE)
> a<-table.element(a,'S.D.',header=TRUE)
> a<-table.element(a,'T-STAT
H0: parameter = 0',header=TRUE)
> a<-table.element(a,'2-tail p-value',header=TRUE)
> a<-table.element(a,'1-tail p-value',header=TRUE)
> a<-table.row.end(a)
> for (i in 1:k){
+ a<-table.row.start(a)
+ a<-table.element(a,rownames(mysum$coefficients)[i],header=TRUE)
+ a<-table.element(a,mysum$coefficients[i,1])
+ a<-table.element(a, round(mysum$coefficients[i,2],6))
+ a<-table.element(a, round(mysum$coefficients[i,3],4))
+ a<-table.element(a, round(mysum$coefficients[i,4],6))
+ a<-table.element(a, round(mysum$coefficients[i,4]/2,6))
+ a<-table.row.end(a)
+ }
> a<-table.end(a)
> table.save(a,file="/var/www/html/rcomp/tmp/129ded1258731433.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Regression Statistics', 2, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple R',1,TRUE)
> a<-table.element(a, sqrt(mysum$r.squared))
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'R-squared',1,TRUE)
> a<-table.element(a, mysum$r.squared)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Adjusted R-squared',1,TRUE)
> a<-table.element(a, mysum$adj.r.squared)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (value)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[1])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (DF numerator)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[2])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (DF denominator)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[3])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'p-value',1,TRUE)
> a<-table.element(a, 1-pf(mysum$fstatistic[1],mysum$fstatistic[2],mysum$fstatistic[3]))
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Residual Statistics', 2, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Residual Standard Deviation',1,TRUE)
> a<-table.element(a, mysum$sigma)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Sum Squared Residuals',1,TRUE)
> a<-table.element(a, sum(myerror*myerror))
> a<-table.row.end(a)
> a<-table.end(a)
> table.save(a,file="/var/www/html/rcomp/tmp/13utii1258731433.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Actuals, Interpolation, and Residuals', 4, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Time or Index', 1, TRUE)
> a<-table.element(a, 'Actuals', 1, TRUE)
> a<-table.element(a, 'Interpolation
Forecast', 1, TRUE)
> a<-table.element(a, 'Residuals
Prediction Error', 1, TRUE)
> a<-table.row.end(a)
> for (i in 1:n) {
+ a<-table.row.start(a)
+ a<-table.element(a,i, 1, TRUE)
+ a<-table.element(a,x[i])
+ a<-table.element(a,x[i]-mysum$resid[i])
+ a<-table.element(a,mysum$resid[i])
+ a<-table.row.end(a)
+ }
> a<-table.end(a)
> table.save(a,file="/var/www/html/rcomp/tmp/14fx761258731433.tab")
> if (n > n25) {
+ a<-table.start()
+ a<-table.row.start(a)
+ a<-table.element(a,'Goldfeld-Quandt test for Heteroskedasticity',4,TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'p-values',header=TRUE)
+ a<-table.element(a,'Alternative Hypothesis',3,header=TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'breakpoint index',header=TRUE)
+ a<-table.element(a,'greater',header=TRUE)
+ a<-table.element(a,'2-sided',header=TRUE)
+ a<-table.element(a,'less',header=TRUE)
+ a<-table.row.end(a)
+ for (mypoint in kp3:nmkm3) {
+ a<-table.row.start(a)
+ a<-table.element(a,mypoint,header=TRUE)
+ a<-table.element(a,gqarr[mypoint-kp3+1,1])
+ a<-table.element(a,gqarr[mypoint-kp3+1,2])
+ a<-table.element(a,gqarr[mypoint-kp3+1,3])
+ a<-table.row.end(a)
+ }
+ a<-table.end(a)
+ table.save(a,file="/var/www/html/rcomp/tmp/15nob71258731433.tab")
+ a<-table.start()
+ a<-table.row.start(a)
+ a<-table.element(a,'Meta Analysis of Goldfeld-Quandt test for Heteroskedasticity',4,TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'Description',header=TRUE)
+ a<-table.element(a,'# significant tests',header=TRUE)
+ a<-table.element(a,'% significant tests',header=TRUE)
+ a<-table.element(a,'OK/NOK',header=TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'1% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant1)
+ a<-table.element(a,numsignificant1/numgqtests)
+ if (numsignificant1/numgqtests < 0.01) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'5% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant5)
+ a<-table.element(a,numsignificant5/numgqtests)
+ if (numsignificant5/numgqtests < 0.05) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'10% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant10)
+ a<-table.element(a,numsignificant10/numgqtests)
+ if (numsignificant10/numgqtests < 0.1) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.end(a)
+ table.save(a,file="/var/www/html/rcomp/tmp/16244e1258731433.tab")
+ }
>
> system("convert tmp/1pwh31258731432.ps tmp/1pwh31258731432.png")
> system("convert tmp/23r791258731432.ps tmp/23r791258731432.png")
> system("convert tmp/3tj1g1258731432.ps tmp/3tj1g1258731432.png")
> system("convert tmp/4vopj1258731432.ps tmp/4vopj1258731432.png")
> system("convert tmp/5fimt1258731432.ps tmp/5fimt1258731432.png")
> system("convert tmp/69yzh1258731432.ps tmp/69yzh1258731432.png")
> system("convert tmp/7kbxq1258731432.ps tmp/7kbxq1258731432.png")
> system("convert tmp/89zjc1258731432.ps tmp/89zjc1258731432.png")
> system("convert tmp/99zx01258731432.ps tmp/99zx01258731432.png")
> system("convert tmp/10hu6a1258731432.ps tmp/10hu6a1258731432.png")
>
>
> proc.time()
user system elapsed
2.361 1.564 2.745