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.
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(493,0,481,0,462,0,457,0,442,0,439,0,488,0,521,0,501,0,485,0,464,0,460,0,467,0,460,0,448,0,443,0,436,0,431,0,484,0,510,0,513,0,503,0,471,0,471,0,476,0,475,0,470,0,461,0,455,0,456,0,517,0,525,0,523,0,519,0,509,0,512,0,519,1,517,1,510,1,509,1,501,1,507,1,569,1,580,1,578,1,565,1,547,1,555,1,562,1,561,1,555,1,544,1,537,1,543,1,594,1,611,1,613,1,611,1,594,1,595,1,591,1,589,1,584,1,573,1,567,1,569,1,621,1,629,1,628,1,612,1,595,1,597,1,593,1,590,1,580,1,574,1,573,1,573,1,620,1,626,1,620,1,588,1,566,1,557,1,561,1,549,1,532,1,526,1,511,1,499,1,555,1,565,1,542,1,527,1,510,1,514,1,517,1,508,1,493,1,490,1,469,1,478,1),dim=c(2,102),dimnames=list(c('Aantal_werklozen_(*1000)','dummyvariabele'),1:102))
> y <- array(NA,dim=c(2,102),dimnames=list(c('Aantal_werklozen_(*1000)','dummyvariabele'),1:102))
> 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
Aantal_werklozen_(*1000) dummyvariabele M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 493 0 1 0 0 0 0 0 0 0 0 0 0
2 481 0 0 1 0 0 0 0 0 0 0 0 0
3 462 0 0 0 1 0 0 0 0 0 0 0 0
4 457 0 0 0 0 1 0 0 0 0 0 0 0
5 442 0 0 0 0 0 1 0 0 0 0 0 0
6 439 0 0 0 0 0 0 1 0 0 0 0 0
7 488 0 0 0 0 0 0 0 1 0 0 0 0
8 521 0 0 0 0 0 0 0 0 1 0 0 0
9 501 0 0 0 0 0 0 0 0 0 1 0 0
10 485 0 0 0 0 0 0 0 0 0 0 1 0
11 464 0 0 0 0 0 0 0 0 0 0 0 1
12 460 0 0 0 0 0 0 0 0 0 0 0 0
13 467 0 1 0 0 0 0 0 0 0 0 0 0
14 460 0 0 1 0 0 0 0 0 0 0 0 0
15 448 0 0 0 1 0 0 0 0 0 0 0 0
16 443 0 0 0 0 1 0 0 0 0 0 0 0
17 436 0 0 0 0 0 1 0 0 0 0 0 0
18 431 0 0 0 0 0 0 1 0 0 0 0 0
19 484 0 0 0 0 0 0 0 1 0 0 0 0
20 510 0 0 0 0 0 0 0 0 1 0 0 0
21 513 0 0 0 0 0 0 0 0 0 1 0 0
22 503 0 0 0 0 0 0 0 0 0 0 1 0
23 471 0 0 0 0 0 0 0 0 0 0 0 1
24 471 0 0 0 0 0 0 0 0 0 0 0 0
25 476 0 1 0 0 0 0 0 0 0 0 0 0
26 475 0 0 1 0 0 0 0 0 0 0 0 0
27 470 0 0 0 1 0 0 0 0 0 0 0 0
28 461 0 0 0 0 1 0 0 0 0 0 0 0
29 455 0 0 0 0 0 1 0 0 0 0 0 0
30 456 0 0 0 0 0 0 1 0 0 0 0 0
31 517 0 0 0 0 0 0 0 1 0 0 0 0
32 525 0 0 0 0 0 0 0 0 1 0 0 0
33 523 0 0 0 0 0 0 0 0 0 1 0 0
34 519 0 0 0 0 0 0 0 0 0 0 1 0
35 509 0 0 0 0 0 0 0 0 0 0 0 1
36 512 0 0 0 0 0 0 0 0 0 0 0 0
37 519 1 1 0 0 0 0 0 0 0 0 0 0
38 517 1 0 1 0 0 0 0 0 0 0 0 0
39 510 1 0 0 1 0 0 0 0 0 0 0 0
40 509 1 0 0 0 1 0 0 0 0 0 0 0
41 501 1 0 0 0 0 1 0 0 0 0 0 0
42 507 1 0 0 0 0 0 1 0 0 0 0 0
43 569 1 0 0 0 0 0 0 1 0 0 0 0
44 580 1 0 0 0 0 0 0 0 1 0 0 0
45 578 1 0 0 0 0 0 0 0 0 1 0 0
46 565 1 0 0 0 0 0 0 0 0 0 1 0
47 547 1 0 0 0 0 0 0 0 0 0 0 1
48 555 1 0 0 0 0 0 0 0 0 0 0 0
49 562 1 1 0 0 0 0 0 0 0 0 0 0
50 561 1 0 1 0 0 0 0 0 0 0 0 0
51 555 1 0 0 1 0 0 0 0 0 0 0 0
52 544 1 0 0 0 1 0 0 0 0 0 0 0
53 537 1 0 0 0 0 1 0 0 0 0 0 0
54 543 1 0 0 0 0 0 1 0 0 0 0 0
55 594 1 0 0 0 0 0 0 1 0 0 0 0
56 611 1 0 0 0 0 0 0 0 1 0 0 0
57 613 1 0 0 0 0 0 0 0 0 1 0 0
58 611 1 0 0 0 0 0 0 0 0 0 1 0
59 594 1 0 0 0 0 0 0 0 0 0 0 1
60 595 1 0 0 0 0 0 0 0 0 0 0 0
61 591 1 1 0 0 0 0 0 0 0 0 0 0
62 589 1 0 1 0 0 0 0 0 0 0 0 0
63 584 1 0 0 1 0 0 0 0 0 0 0 0
64 573 1 0 0 0 1 0 0 0 0 0 0 0
65 567 1 0 0 0 0 1 0 0 0 0 0 0
66 569 1 0 0 0 0 0 1 0 0 0 0 0
67 621 1 0 0 0 0 0 0 1 0 0 0 0
68 629 1 0 0 0 0 0 0 0 1 0 0 0
69 628 1 0 0 0 0 0 0 0 0 1 0 0
70 612 1 0 0 0 0 0 0 0 0 0 1 0
71 595 1 0 0 0 0 0 0 0 0 0 0 1
72 597 1 0 0 0 0 0 0 0 0 0 0 0
73 593 1 1 0 0 0 0 0 0 0 0 0 0
74 590 1 0 1 0 0 0 0 0 0 0 0 0
75 580 1 0 0 1 0 0 0 0 0 0 0 0
76 574 1 0 0 0 1 0 0 0 0 0 0 0
77 573 1 0 0 0 0 1 0 0 0 0 0 0
78 573 1 0 0 0 0 0 1 0 0 0 0 0
79 620 1 0 0 0 0 0 0 1 0 0 0 0
80 626 1 0 0 0 0 0 0 0 1 0 0 0
81 620 1 0 0 0 0 0 0 0 0 1 0 0
82 588 1 0 0 0 0 0 0 0 0 0 1 0
83 566 1 0 0 0 0 0 0 0 0 0 0 1
84 557 1 0 0 0 0 0 0 0 0 0 0 0
85 561 1 1 0 0 0 0 0 0 0 0 0 0
86 549 1 0 1 0 0 0 0 0 0 0 0 0
87 532 1 0 0 1 0 0 0 0 0 0 0 0
88 526 1 0 0 0 1 0 0 0 0 0 0 0
89 511 1 0 0 0 0 1 0 0 0 0 0 0
90 499 1 0 0 0 0 0 1 0 0 0 0 0
91 555 1 0 0 0 0 0 0 1 0 0 0 0
92 565 1 0 0 0 0 0 0 0 1 0 0 0
93 542 1 0 0 0 0 0 0 0 0 1 0 0
94 527 1 0 0 0 0 0 0 0 0 0 1 0
95 510 1 0 0 0 0 0 0 0 0 0 0 1
96 514 1 0 0 0 0 0 0 0 0 0 0 0
97 517 1 1 0 0 0 0 0 0 0 0 0 0
98 508 1 0 1 0 0 0 0 0 0 0 0 0
99 493 1 0 0 1 0 0 0 0 0 0 0 0
100 490 1 0 0 0 1 0 0 0 0 0 0 0
101 469 1 0 0 0 0 1 0 0 0 0 0 0
102 478 1 0 0 0 0 0 1 0 0 0 0 0
t
1 1
2 2
3 3
4 4
5 5
6 6
7 7
8 8
9 9
10 10
11 11
12 12
13 13
14 14
15 15
16 16
17 17
18 18
19 19
20 20
21 21
22 22
23 23
24 24
25 25
26 26
27 27
28 28
29 29
30 30
31 31
32 32
33 33
34 34
35 35
36 36
37 37
38 38
39 39
40 40
41 41
42 42
43 43
44 44
45 45
46 46
47 47
48 48
49 49
50 50
51 51
52 52
53 53
54 54
55 55
56 56
57 57
58 58
59 59
60 60
61 61
62 62
63 63
64 64
65 65
66 66
67 67
68 68
69 69
70 70
71 71
72 72
73 73
74 74
75 75
76 76
77 77
78 78
79 79
80 80
81 81
82 82
83 83
84 84
85 85
86 86
87 87
88 88
89 89
90 90
91 91
92 92
93 93
94 94
95 95
96 96
97 97
98 98
99 99
100 100
101 101
102 102
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) dummyvariabele M1 M2 M3
487.0749 97.6558 -7.1278 -12.2855 -22.6654
M4 M5 M6 M7 M8
-28.7119 -37.9807 -37.2495 21.9412 37.1030
M9 M10 M11 t
31.2647 18.0515 -0.9118 -0.2868
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-48.827 -19.428 0.949 24.543 48.330
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 487.0749 11.3310 42.986 < 2e-16 ***
dummyvariabele 97.6558 10.5460 9.260 1.20e-14 ***
M1 -7.1278 13.8401 -0.515 0.60784
M2 -12.2855 13.8260 -0.889 0.37665
M3 -22.6654 13.8141 -1.641 0.10442
M4 -28.7119 13.8043 -2.080 0.04044 *
M5 -37.9807 13.7965 -2.753 0.00717 **
M6 -37.2495 13.7910 -2.701 0.00829 **
M7 21.9412 14.2094 1.544 0.12614
M8 37.1030 14.2001 2.613 0.01056 *
M9 31.2647 14.1929 2.203 0.03022 *
M10 18.0515 14.1877 1.272 0.20661
M11 -0.9118 14.1846 -0.064 0.94889
t -0.2868 0.1713 -1.674 0.09767 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 28.37 on 88 degrees of freedom
Multiple R-squared: 0.7547, Adjusted R-squared: 0.7185
F-statistic: 20.83 on 13 and 88 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,] 1.125154e-02 2.250308e-02 0.98874846
[2,] 2.836508e-03 5.673016e-03 0.99716349
[3,] 9.987605e-04 1.997521e-03 0.99900124
[4,] 1.837770e-04 3.675540e-04 0.99981622
[5,] 5.816958e-04 1.163392e-03 0.99941830
[6,] 1.114090e-03 2.228180e-03 0.99888591
[7,] 5.246224e-04 1.049245e-03 0.99947538
[8,] 2.979160e-04 5.958321e-04 0.99970208
[9,] 9.211393e-05 1.842279e-04 0.99990789
[10,] 3.682364e-05 7.364728e-05 0.99996318
[11,] 2.946109e-05 5.892217e-05 0.99997054
[12,] 1.342347e-05 2.684694e-05 0.99998658
[13,] 7.828869e-06 1.565774e-05 0.99999217
[14,] 5.969996e-06 1.193999e-05 0.99999403
[15,] 9.701594e-06 1.940319e-05 0.99999030
[16,] 3.340092e-06 6.680183e-06 0.99999666
[17,] 1.336952e-06 2.673904e-06 0.99999866
[18,] 8.520353e-07 1.704071e-06 0.99999915
[19,] 2.282823e-06 4.565646e-06 0.99999772
[20,] 6.010060e-06 1.202012e-05 0.99999399
[21,] 3.690308e-06 7.380616e-06 0.99999631
[22,] 2.403913e-06 4.807827e-06 0.99999760
[23,] 1.715236e-06 3.430472e-06 0.99999828
[24,] 1.320475e-06 2.640950e-06 0.99999868
[25,] 1.084867e-06 2.169734e-06 0.99999892
[26,] 1.201028e-06 2.402056e-06 0.99999880
[27,] 1.889632e-06 3.779263e-06 0.99999811
[28,] 1.852413e-06 3.704826e-06 0.99999815
[29,] 2.051055e-06 4.102109e-06 0.99999795
[30,] 2.219801e-06 4.439602e-06 0.99999778
[31,] 3.130767e-06 6.261534e-06 0.99999687
[32,] 5.551844e-06 1.110369e-05 0.99999445
[33,] 1.226296e-05 2.452592e-05 0.99998774
[34,] 3.153695e-05 6.307391e-05 0.99996846
[35,] 9.134616e-05 1.826923e-04 0.99990865
[36,] 2.334291e-04 4.668582e-04 0.99976657
[37,] 6.821089e-04 1.364218e-03 0.99931789
[38,] 2.428755e-03 4.857511e-03 0.99757124
[39,] 9.741666e-03 1.948333e-02 0.99025833
[40,] 2.865698e-02 5.731396e-02 0.97134302
[41,] 6.636441e-02 1.327288e-01 0.93363559
[42,] 1.041431e-01 2.082862e-01 0.89585692
[43,] 1.610339e-01 3.220679e-01 0.83896607
[44,] 2.229280e-01 4.458561e-01 0.77707197
[45,] 3.077103e-01 6.154207e-01 0.69228967
[46,] 3.909988e-01 7.819977e-01 0.60900116
[47,] 4.462591e-01 8.925183e-01 0.55374087
[48,] 5.626401e-01 8.747197e-01 0.43735986
[49,] 6.613906e-01 6.772187e-01 0.33860937
[50,] 7.644193e-01 4.711614e-01 0.23558068
[51,] 8.468082e-01 3.063835e-01 0.15319176
[52,] 9.257552e-01 1.484896e-01 0.07424478
[53,] 9.375289e-01 1.249421e-01 0.06247107
[54,] 9.301823e-01 1.396355e-01 0.06981774
[55,] 9.179966e-01 1.640068e-01 0.08200341
[56,] 8.933095e-01 2.133810e-01 0.10669048
[57,] 9.040984e-01 1.918033e-01 0.09590163
[58,] 8.880796e-01 2.238409e-01 0.11192043
[59,] 8.509943e-01 2.980113e-01 0.14900567
[60,] 8.159814e-01 3.680372e-01 0.18401859
[61,] 7.574940e-01 4.850120e-01 0.24250600
[62,] 7.048826e-01 5.902348e-01 0.29511741
[63,] 6.875663e-01 6.248673e-01 0.31243367
[64,] 6.578429e-01 6.843142e-01 0.34215709
[65,] 8.984936e-01 2.030129e-01 0.10150644
[66,] 9.483416e-01 1.033168e-01 0.05165840
[67,] 9.721981e-01 5.560388e-02 0.02780194
[68,] 9.495788e-01 1.008425e-01 0.05042124
[69,] 9.074655e-01 1.850690e-01 0.09253450
> postscript(file="/var/www/html/freestat/rcomp/tmp/13q8q1228658906.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/freestat/rcomp/tmp/2s2lm1228658906.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/freestat/rcomp/tmp/3xb3n1228658906.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/freestat/rcomp/tmp/4j98n1228658906.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/freestat/rcomp/tmp/54o001228658906.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 = 102
Frequency = 1
1 2 3 4 5 6
13.3396567 6.7841012 -1.5492322 -0.2158988 -5.6603433 -9.1047877
7 8 9 10 11 12
-19.0088076 -0.8838076 -14.7588076 -17.2588076 -19.0088076 -23.6338076
13 14 15 16 17 18
-9.2192864 -10.7748419 -12.1081752 -10.7748419 -8.2192864 -13.6637308
19 20 21 22 23 24
-19.5677507 -8.4427507 0.6822493 4.1822493 -8.5677507 -9.1927507
25 26 27 28 29 30
3.2217706 7.6662150 13.3328817 10.6662150 14.2217706 14.7773261
31 32 33 34 35 36
16.8733062 9.9983062 14.1233062 23.6233062 32.8733062 35.2483062
37 38 39 40 41 42
-47.9929991 -44.5485547 -40.8818880 -35.5485547 -33.9929991 -28.4374435
43 44 45 46 47 48
-25.3414634 -29.2164634 -25.0914634 -24.5914634 -23.3414634 -15.9664634
49 50 51 52 53 54
-1.5519422 2.8925023 7.5591689 2.8925023 5.4480578 11.0036134
55 56 57 58 59 60
3.0995935 5.2245935 13.3495935 24.8495935 27.0995935 27.4745935
61 62 63 64 65 66
30.8891147 34.3335592 40.0002258 35.3335592 38.8891147 40.4446703
67 68 69 70 71 72
33.5406504 26.6656504 31.7906504 29.2906504 31.5406504 32.9156504
73 74 75 76 77 78
36.3301716 38.7746161 39.4412827 39.7746161 48.3301716 47.8857272
79 80 81 82 83 84
35.9817073 27.1067073 27.2317073 8.7317073 5.9817073 -3.6432927
85 86 87 88 89 90
7.7712285 1.2156730 -5.1176603 -4.7843270 -10.2287715 -22.6732159
91 92 93 94 95 96
-25.5772358 -30.4522358 -47.3272358 -48.8272358 -46.5772358 -43.2022358
97 98 99 100 101 102
-32.7877145 -36.3432701 -40.6766034 -37.3432701 -48.7877145 -40.2321590
> postscript(file="/var/www/html/freestat/rcomp/tmp/6pwc21228658906.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 = 102
Frequency = 1
lag(myerror, k = 1) myerror
0 13.3396567 NA
1 6.7841012 13.3396567
2 -1.5492322 6.7841012
3 -0.2158988 -1.5492322
4 -5.6603433 -0.2158988
5 -9.1047877 -5.6603433
6 -19.0088076 -9.1047877
7 -0.8838076 -19.0088076
8 -14.7588076 -0.8838076
9 -17.2588076 -14.7588076
10 -19.0088076 -17.2588076
11 -23.6338076 -19.0088076
12 -9.2192864 -23.6338076
13 -10.7748419 -9.2192864
14 -12.1081752 -10.7748419
15 -10.7748419 -12.1081752
16 -8.2192864 -10.7748419
17 -13.6637308 -8.2192864
18 -19.5677507 -13.6637308
19 -8.4427507 -19.5677507
20 0.6822493 -8.4427507
21 4.1822493 0.6822493
22 -8.5677507 4.1822493
23 -9.1927507 -8.5677507
24 3.2217706 -9.1927507
25 7.6662150 3.2217706
26 13.3328817 7.6662150
27 10.6662150 13.3328817
28 14.2217706 10.6662150
29 14.7773261 14.2217706
30 16.8733062 14.7773261
31 9.9983062 16.8733062
32 14.1233062 9.9983062
33 23.6233062 14.1233062
34 32.8733062 23.6233062
35 35.2483062 32.8733062
36 -47.9929991 35.2483062
37 -44.5485547 -47.9929991
38 -40.8818880 -44.5485547
39 -35.5485547 -40.8818880
40 -33.9929991 -35.5485547
41 -28.4374435 -33.9929991
42 -25.3414634 -28.4374435
43 -29.2164634 -25.3414634
44 -25.0914634 -29.2164634
45 -24.5914634 -25.0914634
46 -23.3414634 -24.5914634
47 -15.9664634 -23.3414634
48 -1.5519422 -15.9664634
49 2.8925023 -1.5519422
50 7.5591689 2.8925023
51 2.8925023 7.5591689
52 5.4480578 2.8925023
53 11.0036134 5.4480578
54 3.0995935 11.0036134
55 5.2245935 3.0995935
56 13.3495935 5.2245935
57 24.8495935 13.3495935
58 27.0995935 24.8495935
59 27.4745935 27.0995935
60 30.8891147 27.4745935
61 34.3335592 30.8891147
62 40.0002258 34.3335592
63 35.3335592 40.0002258
64 38.8891147 35.3335592
65 40.4446703 38.8891147
66 33.5406504 40.4446703
67 26.6656504 33.5406504
68 31.7906504 26.6656504
69 29.2906504 31.7906504
70 31.5406504 29.2906504
71 32.9156504 31.5406504
72 36.3301716 32.9156504
73 38.7746161 36.3301716
74 39.4412827 38.7746161
75 39.7746161 39.4412827
76 48.3301716 39.7746161
77 47.8857272 48.3301716
78 35.9817073 47.8857272
79 27.1067073 35.9817073
80 27.2317073 27.1067073
81 8.7317073 27.2317073
82 5.9817073 8.7317073
83 -3.6432927 5.9817073
84 7.7712285 -3.6432927
85 1.2156730 7.7712285
86 -5.1176603 1.2156730
87 -4.7843270 -5.1176603
88 -10.2287715 -4.7843270
89 -22.6732159 -10.2287715
90 -25.5772358 -22.6732159
91 -30.4522358 -25.5772358
92 -47.3272358 -30.4522358
93 -48.8272358 -47.3272358
94 -46.5772358 -48.8272358
95 -43.2022358 -46.5772358
96 -32.7877145 -43.2022358
97 -36.3432701 -32.7877145
98 -40.6766034 -36.3432701
99 -37.3432701 -40.6766034
100 -48.7877145 -37.3432701
101 -40.2321590 -48.7877145
102 NA -40.2321590
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 6.7841012 13.3396567
[2,] -1.5492322 6.7841012
[3,] -0.2158988 -1.5492322
[4,] -5.6603433 -0.2158988
[5,] -9.1047877 -5.6603433
[6,] -19.0088076 -9.1047877
[7,] -0.8838076 -19.0088076
[8,] -14.7588076 -0.8838076
[9,] -17.2588076 -14.7588076
[10,] -19.0088076 -17.2588076
[11,] -23.6338076 -19.0088076
[12,] -9.2192864 -23.6338076
[13,] -10.7748419 -9.2192864
[14,] -12.1081752 -10.7748419
[15,] -10.7748419 -12.1081752
[16,] -8.2192864 -10.7748419
[17,] -13.6637308 -8.2192864
[18,] -19.5677507 -13.6637308
[19,] -8.4427507 -19.5677507
[20,] 0.6822493 -8.4427507
[21,] 4.1822493 0.6822493
[22,] -8.5677507 4.1822493
[23,] -9.1927507 -8.5677507
[24,] 3.2217706 -9.1927507
[25,] 7.6662150 3.2217706
[26,] 13.3328817 7.6662150
[27,] 10.6662150 13.3328817
[28,] 14.2217706 10.6662150
[29,] 14.7773261 14.2217706
[30,] 16.8733062 14.7773261
[31,] 9.9983062 16.8733062
[32,] 14.1233062 9.9983062
[33,] 23.6233062 14.1233062
[34,] 32.8733062 23.6233062
[35,] 35.2483062 32.8733062
[36,] -47.9929991 35.2483062
[37,] -44.5485547 -47.9929991
[38,] -40.8818880 -44.5485547
[39,] -35.5485547 -40.8818880
[40,] -33.9929991 -35.5485547
[41,] -28.4374435 -33.9929991
[42,] -25.3414634 -28.4374435
[43,] -29.2164634 -25.3414634
[44,] -25.0914634 -29.2164634
[45,] -24.5914634 -25.0914634
[46,] -23.3414634 -24.5914634
[47,] -15.9664634 -23.3414634
[48,] -1.5519422 -15.9664634
[49,] 2.8925023 -1.5519422
[50,] 7.5591689 2.8925023
[51,] 2.8925023 7.5591689
[52,] 5.4480578 2.8925023
[53,] 11.0036134 5.4480578
[54,] 3.0995935 11.0036134
[55,] 5.2245935 3.0995935
[56,] 13.3495935 5.2245935
[57,] 24.8495935 13.3495935
[58,] 27.0995935 24.8495935
[59,] 27.4745935 27.0995935
[60,] 30.8891147 27.4745935
[61,] 34.3335592 30.8891147
[62,] 40.0002258 34.3335592
[63,] 35.3335592 40.0002258
[64,] 38.8891147 35.3335592
[65,] 40.4446703 38.8891147
[66,] 33.5406504 40.4446703
[67,] 26.6656504 33.5406504
[68,] 31.7906504 26.6656504
[69,] 29.2906504 31.7906504
[70,] 31.5406504 29.2906504
[71,] 32.9156504 31.5406504
[72,] 36.3301716 32.9156504
[73,] 38.7746161 36.3301716
[74,] 39.4412827 38.7746161
[75,] 39.7746161 39.4412827
[76,] 48.3301716 39.7746161
[77,] 47.8857272 48.3301716
[78,] 35.9817073 47.8857272
[79,] 27.1067073 35.9817073
[80,] 27.2317073 27.1067073
[81,] 8.7317073 27.2317073
[82,] 5.9817073 8.7317073
[83,] -3.6432927 5.9817073
[84,] 7.7712285 -3.6432927
[85,] 1.2156730 7.7712285
[86,] -5.1176603 1.2156730
[87,] -4.7843270 -5.1176603
[88,] -10.2287715 -4.7843270
[89,] -22.6732159 -10.2287715
[90,] -25.5772358 -22.6732159
[91,] -30.4522358 -25.5772358
[92,] -47.3272358 -30.4522358
[93,] -48.8272358 -47.3272358
[94,] -46.5772358 -48.8272358
[95,] -43.2022358 -46.5772358
[96,] -32.7877145 -43.2022358
[97,] -36.3432701 -32.7877145
[98,] -40.6766034 -36.3432701
[99,] -37.3432701 -40.6766034
[100,] -48.7877145 -37.3432701
[101,] -40.2321590 -48.7877145
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 6.7841012 13.3396567
2 -1.5492322 6.7841012
3 -0.2158988 -1.5492322
4 -5.6603433 -0.2158988
5 -9.1047877 -5.6603433
6 -19.0088076 -9.1047877
7 -0.8838076 -19.0088076
8 -14.7588076 -0.8838076
9 -17.2588076 -14.7588076
10 -19.0088076 -17.2588076
11 -23.6338076 -19.0088076
12 -9.2192864 -23.6338076
13 -10.7748419 -9.2192864
14 -12.1081752 -10.7748419
15 -10.7748419 -12.1081752
16 -8.2192864 -10.7748419
17 -13.6637308 -8.2192864
18 -19.5677507 -13.6637308
19 -8.4427507 -19.5677507
20 0.6822493 -8.4427507
21 4.1822493 0.6822493
22 -8.5677507 4.1822493
23 -9.1927507 -8.5677507
24 3.2217706 -9.1927507
25 7.6662150 3.2217706
26 13.3328817 7.6662150
27 10.6662150 13.3328817
28 14.2217706 10.6662150
29 14.7773261 14.2217706
30 16.8733062 14.7773261
31 9.9983062 16.8733062
32 14.1233062 9.9983062
33 23.6233062 14.1233062
34 32.8733062 23.6233062
35 35.2483062 32.8733062
36 -47.9929991 35.2483062
37 -44.5485547 -47.9929991
38 -40.8818880 -44.5485547
39 -35.5485547 -40.8818880
40 -33.9929991 -35.5485547
41 -28.4374435 -33.9929991
42 -25.3414634 -28.4374435
43 -29.2164634 -25.3414634
44 -25.0914634 -29.2164634
45 -24.5914634 -25.0914634
46 -23.3414634 -24.5914634
47 -15.9664634 -23.3414634
48 -1.5519422 -15.9664634
49 2.8925023 -1.5519422
50 7.5591689 2.8925023
51 2.8925023 7.5591689
52 5.4480578 2.8925023
53 11.0036134 5.4480578
54 3.0995935 11.0036134
55 5.2245935 3.0995935
56 13.3495935 5.2245935
57 24.8495935 13.3495935
58 27.0995935 24.8495935
59 27.4745935 27.0995935
60 30.8891147 27.4745935
61 34.3335592 30.8891147
62 40.0002258 34.3335592
63 35.3335592 40.0002258
64 38.8891147 35.3335592
65 40.4446703 38.8891147
66 33.5406504 40.4446703
67 26.6656504 33.5406504
68 31.7906504 26.6656504
69 29.2906504 31.7906504
70 31.5406504 29.2906504
71 32.9156504 31.5406504
72 36.3301716 32.9156504
73 38.7746161 36.3301716
74 39.4412827 38.7746161
75 39.7746161 39.4412827
76 48.3301716 39.7746161
77 47.8857272 48.3301716
78 35.9817073 47.8857272
79 27.1067073 35.9817073
80 27.2317073 27.1067073
81 8.7317073 27.2317073
82 5.9817073 8.7317073
83 -3.6432927 5.9817073
84 7.7712285 -3.6432927
85 1.2156730 7.7712285
86 -5.1176603 1.2156730
87 -4.7843270 -5.1176603
88 -10.2287715 -4.7843270
89 -22.6732159 -10.2287715
90 -25.5772358 -22.6732159
91 -30.4522358 -25.5772358
92 -47.3272358 -30.4522358
93 -48.8272358 -47.3272358
94 -46.5772358 -48.8272358
95 -43.2022358 -46.5772358
96 -32.7877145 -43.2022358
97 -36.3432701 -32.7877145
98 -40.6766034 -36.3432701
99 -37.3432701 -40.6766034
100 -48.7877145 -37.3432701
101 -40.2321590 -48.7877145
> 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/76vmx1228658906.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/freestat/rcomp/tmp/87dio1228658906.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/freestat/rcomp/tmp/9ndfa1228658906.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/freestat/rcomp/tmp/10wf8m1228658906.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/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/1128ol1228658906.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/12ryly1228658907.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/13lkr21228658907.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/14nt0h1228658907.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/15hj991228658907.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/16kz3h1228658907.tab")
+ }
>
> system("convert tmp/13q8q1228658906.ps tmp/13q8q1228658906.png")
> system("convert tmp/2s2lm1228658906.ps tmp/2s2lm1228658906.png")
> system("convert tmp/3xb3n1228658906.ps tmp/3xb3n1228658906.png")
> system("convert tmp/4j98n1228658906.ps tmp/4j98n1228658906.png")
> system("convert tmp/54o001228658906.ps tmp/54o001228658906.png")
> system("convert tmp/6pwc21228658906.ps tmp/6pwc21228658906.png")
> system("convert tmp/76vmx1228658906.ps tmp/76vmx1228658906.png")
> system("convert tmp/87dio1228658906.ps tmp/87dio1228658906.png")
> system("convert tmp/9ndfa1228658906.ps tmp/9ndfa1228658906.png")
> system("convert tmp/10wf8m1228658906.ps tmp/10wf8m1228658906.png")
>
>
> proc.time()
user system elapsed
4.357 2.530 4.691