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(524,0,552,0,532,0,511,0,492,0,492,0,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,0,517,0,510,0,509,0,501,0,507,0,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),dim=c(2,96),dimnames=list(c('Aantal_werklozen_(*1000)','Xt'),1:96))
> y <- array(NA,dim=c(2,96),dimnames=list(c('Aantal_werklozen_(*1000)','Xt'),1:96))
> 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) Xt M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 524 0 1 0 0 0 0 0 0 0 0 0 0 1
2 552 0 0 1 0 0 0 0 0 0 0 0 0 2
3 532 0 0 0 1 0 0 0 0 0 0 0 0 3
4 511 0 0 0 0 1 0 0 0 0 0 0 0 4
5 492 0 0 0 0 0 1 0 0 0 0 0 0 5
6 492 0 0 0 0 0 0 1 0 0 0 0 0 6
7 493 0 0 0 0 0 0 0 1 0 0 0 0 7
8 481 0 0 0 0 0 0 0 0 1 0 0 0 8
9 462 0 0 0 0 0 0 0 0 0 1 0 0 9
10 457 0 0 0 0 0 0 0 0 0 0 1 0 10
11 442 0 0 0 0 0 0 0 0 0 0 0 1 11
12 439 0 0 0 0 0 0 0 0 0 0 0 0 12
13 488 0 1 0 0 0 0 0 0 0 0 0 0 13
14 521 0 0 1 0 0 0 0 0 0 0 0 0 14
15 501 0 0 0 1 0 0 0 0 0 0 0 0 15
16 485 0 0 0 0 1 0 0 0 0 0 0 0 16
17 464 0 0 0 0 0 1 0 0 0 0 0 0 17
18 460 0 0 0 0 0 0 1 0 0 0 0 0 18
19 467 0 0 0 0 0 0 0 1 0 0 0 0 19
20 460 0 0 0 0 0 0 0 0 1 0 0 0 20
21 448 0 0 0 0 0 0 0 0 0 1 0 0 21
22 443 0 0 0 0 0 0 0 0 0 0 1 0 22
23 436 0 0 0 0 0 0 0 0 0 0 0 1 23
24 431 0 0 0 0 0 0 0 0 0 0 0 0 24
25 484 0 1 0 0 0 0 0 0 0 0 0 0 25
26 510 0 0 1 0 0 0 0 0 0 0 0 0 26
27 513 0 0 0 1 0 0 0 0 0 0 0 0 27
28 503 0 0 0 0 1 0 0 0 0 0 0 0 28
29 471 0 0 0 0 0 1 0 0 0 0 0 0 29
30 471 0 0 0 0 0 0 1 0 0 0 0 0 30
31 476 0 0 0 0 0 0 0 1 0 0 0 0 31
32 475 0 0 0 0 0 0 0 0 1 0 0 0 32
33 470 0 0 0 0 0 0 0 0 0 1 0 0 33
34 461 0 0 0 0 0 0 0 0 0 0 1 0 34
35 455 0 0 0 0 0 0 0 0 0 0 0 1 35
36 456 0 0 0 0 0 0 0 0 0 0 0 0 36
37 517 0 1 0 0 0 0 0 0 0 0 0 0 37
38 525 0 0 1 0 0 0 0 0 0 0 0 0 38
39 523 0 0 0 1 0 0 0 0 0 0 0 0 39
40 519 0 0 0 0 1 0 0 0 0 0 0 0 40
41 509 0 0 0 0 0 1 0 0 0 0 0 0 41
42 512 0 0 0 0 0 0 1 0 0 0 0 0 42
43 519 0 0 0 0 0 0 0 1 0 0 0 0 43
44 517 0 0 0 0 0 0 0 0 1 0 0 0 44
45 510 0 0 0 0 0 0 0 0 0 1 0 0 45
46 509 0 0 0 0 0 0 0 0 0 0 1 0 46
47 501 0 0 0 0 0 0 0 0 0 0 0 1 47
48 507 0 0 0 0 0 0 0 0 0 0 0 0 48
49 569 1 1 0 0 0 0 0 0 0 0 0 0 49
50 580 1 0 1 0 0 0 0 0 0 0 0 0 50
51 578 1 0 0 1 0 0 0 0 0 0 0 0 51
52 565 1 0 0 0 1 0 0 0 0 0 0 0 52
53 547 1 0 0 0 0 1 0 0 0 0 0 0 53
54 555 1 0 0 0 0 0 1 0 0 0 0 0 54
55 562 1 0 0 0 0 0 0 1 0 0 0 0 55
56 561 1 0 0 0 0 0 0 0 1 0 0 0 56
57 555 1 0 0 0 0 0 0 0 0 1 0 0 57
58 544 1 0 0 0 0 0 0 0 0 0 1 0 58
59 537 1 0 0 0 0 0 0 0 0 0 0 1 59
60 543 1 0 0 0 0 0 0 0 0 0 0 0 60
61 594 1 1 0 0 0 0 0 0 0 0 0 0 61
62 611 1 0 1 0 0 0 0 0 0 0 0 0 62
63 613 1 0 0 1 0 0 0 0 0 0 0 0 63
64 611 1 0 0 0 1 0 0 0 0 0 0 0 64
65 594 1 0 0 0 0 1 0 0 0 0 0 0 65
66 595 1 0 0 0 0 0 1 0 0 0 0 0 66
67 591 1 0 0 0 0 0 0 1 0 0 0 0 67
68 589 1 0 0 0 0 0 0 0 1 0 0 0 68
69 584 1 0 0 0 0 0 0 0 0 1 0 0 69
70 573 1 0 0 0 0 0 0 0 0 0 1 0 70
71 567 1 0 0 0 0 0 0 0 0 0 0 1 71
72 569 1 0 0 0 0 0 0 0 0 0 0 0 72
73 621 1 1 0 0 0 0 0 0 0 0 0 0 73
74 629 1 0 1 0 0 0 0 0 0 0 0 0 74
75 628 1 0 0 1 0 0 0 0 0 0 0 0 75
76 612 1 0 0 0 1 0 0 0 0 0 0 0 76
77 595 1 0 0 0 0 1 0 0 0 0 0 0 77
78 597 1 0 0 0 0 0 1 0 0 0 0 0 78
79 593 1 0 0 0 0 0 0 1 0 0 0 0 79
80 590 1 0 0 0 0 0 0 0 1 0 0 0 80
81 580 1 0 0 0 0 0 0 0 0 1 0 0 81
82 574 1 0 0 0 0 0 0 0 0 0 1 0 82
83 573 1 0 0 0 0 0 0 0 0 0 0 1 83
84 573 1 0 0 0 0 0 0 0 0 0 0 0 84
85 620 1 1 0 0 0 0 0 0 0 0 0 0 85
86 626 1 0 1 0 0 0 0 0 0 0 0 0 86
87 620 1 0 0 1 0 0 0 0 0 0 0 0 87
88 588 1 0 0 0 1 0 0 0 0 0 0 0 88
89 566 1 0 0 0 0 1 0 0 0 0 0 0 89
90 557 1 0 0 0 0 0 1 0 0 0 0 0 90
91 561 1 0 0 0 0 0 0 1 0 0 0 0 91
92 549 1 0 0 0 0 0 0 0 1 0 0 0 92
93 532 1 0 0 0 0 0 0 0 0 1 0 0 93
94 526 1 0 0 0 0 0 0 0 0 0 1 0 94
95 511 1 0 0 0 0 0 0 0 0 0 0 1 95
96 499 1 0 0 0 0 0 0 0 0 0 0 0 96
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Xt M1 M2 M3 M4
444.0208 68.7708 54.8316 71.5174 65.3281 50.6389
M5 M6 M7 M8 M9 M10
30.6997 30.3854 32.8212 27.3819 16.8177 9.6285
M11 t
1.0642 0.4392
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-55.9583 -16.5833 -0.8646 16.4375 41.8958
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 444.0208 9.3445 47.517 < 2e-16 ***
Xt 68.7708 9.0277 7.618 3.97e-11 ***
M1 54.8316 10.9402 5.012 3.05e-06 ***
M2 71.5174 10.9143 6.553 4.66e-09 ***
M3 65.3281 10.8908 5.998 5.16e-08 ***
M4 50.6389 10.8698 4.659 1.21e-05 ***
M5 30.6997 10.8512 2.829 0.00587 **
M6 30.3854 10.8350 2.804 0.00629 **
M7 32.8212 10.8213 3.033 0.00324 **
M8 27.3819 10.8101 2.533 0.01321 *
M9 16.8177 10.8014 1.557 0.12332
M10 9.6285 10.7951 0.892 0.37504
M11 1.0642 10.7914 0.099 0.92168
t 0.4392 0.1642 2.676 0.00901 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 21.58 on 82 degrees of freedom
Multiple R-squared: 0.8633, Adjusted R-squared: 0.8416
F-statistic: 39.84 on 13 and 82 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,] 3.560278e-03 0.0071205552 0.9964397
[2,] 4.306296e-04 0.0008612593 0.9995694
[3,] 8.736665e-05 0.0001747333 0.9999126
[4,] 6.493674e-05 0.0001298735 0.9999351
[5,] 1.646335e-04 0.0003292670 0.9998354
[6,] 1.425271e-04 0.0002850543 0.9998575
[7,] 3.283937e-04 0.0006567874 0.9996716
[8,] 3.108021e-04 0.0006216042 0.9996892
[9,] 2.095226e-04 0.0004190451 0.9997905
[10,] 7.731907e-05 0.0001546381 0.9999227
[11,] 4.276103e-04 0.0008552206 0.9995724
[12,] 1.976053e-03 0.0039521066 0.9980239
[13,] 1.626588e-03 0.0032531755 0.9983734
[14,] 1.409376e-03 0.0028187523 0.9985906
[15,] 1.175996e-03 0.0023519923 0.9988240
[16,] 1.542859e-03 0.0030857185 0.9984571
[17,] 3.562940e-03 0.0071258800 0.9964371
[18,] 4.707961e-03 0.0094159214 0.9952920
[19,] 7.185943e-03 0.0143718866 0.9928141
[20,] 1.229357e-02 0.0245871318 0.9877064
[21,] 1.704869e-02 0.0340973840 0.9829513
[22,] 1.242892e-02 0.0248578435 0.9875711
[23,] 1.118049e-02 0.0223609815 0.9888195
[24,] 1.280287e-02 0.0256057315 0.9871971
[25,] 2.282122e-02 0.0456424393 0.9771788
[26,] 3.764760e-02 0.0752951932 0.9623524
[27,] 5.398107e-02 0.1079621490 0.9460189
[28,] 7.568829e-02 0.1513765846 0.9243117
[29,] 1.043933e-01 0.2087865675 0.8956067
[30,] 1.404247e-01 0.2808494164 0.8595753
[31,] 1.715203e-01 0.3430406406 0.8284797
[32,] 2.200123e-01 0.4400246866 0.7799877
[33,] 2.179761e-01 0.4359522599 0.7820239
[34,] 2.266927e-01 0.4533854871 0.7733073
[35,] 2.452357e-01 0.4904713133 0.7547643
[36,] 2.666844e-01 0.5333687044 0.7333156
[37,] 3.060665e-01 0.6121329387 0.6939335
[38,] 3.269904e-01 0.6539807964 0.6730096
[39,] 3.287268e-01 0.6574536993 0.6712732
[40,] 3.273895e-01 0.6547789290 0.6726105
[41,] 3.234028e-01 0.6468055602 0.6765972
[42,] 3.379575e-01 0.6759149090 0.6620425
[43,] 3.761240e-01 0.7522479929 0.6238760
[44,] 3.960657e-01 0.7921313917 0.6039343
[45,] 5.467921e-01 0.9064158851 0.4532079
[46,] 6.583830e-01 0.6832339406 0.3416170
[47,] 7.684159e-01 0.4631682792 0.2315841
[48,] 7.727326e-01 0.4545347953 0.2272674
[49,] 7.746933e-01 0.4506133648 0.2253067
[50,] 7.604532e-01 0.4790935132 0.2395468
[51,] 7.569275e-01 0.4861450524 0.2430725
[52,] 7.391753e-01 0.5216494029 0.2608247
[53,] 6.965036e-01 0.6069927632 0.3034964
[54,] 6.756416e-01 0.6487168392 0.3243584
[55,] 6.639275e-01 0.6721449847 0.3360725
[56,] 6.246213e-01 0.7507573940 0.3753787
[57,] 6.973490e-01 0.6053019962 0.3026510
[58,] 7.897044e-01 0.4205911382 0.2102956
[59,] 8.819353e-01 0.2361293168 0.1180647
[60,] 8.814057e-01 0.2371886474 0.1185943
[61,] 8.678756e-01 0.2642488405 0.1321244
[62,] 7.866907e-01 0.4266185277 0.2133093
[63,] 7.803410e-01 0.4393180947 0.2196590
> postscript(file="/var/www/html/rcomp/tmp/15muw1229456596.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/26jle1229456596.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/399np1229456596.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/45kyj1229456596.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/5uaxl1229456596.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 = 96
Frequency = 1
1 2 3 4 5 6
24.7083333 35.5833333 21.3333333 14.5833333 15.0833333 14.9583333
7 8 9 10 11 12
13.0833333 6.0833333 -2.7916667 -1.0416667 -7.9166667 -10.2916667
13 14 15 16 17 18
-16.5625000 -0.6875000 -14.9375000 -16.6875000 -18.1875000 -22.3125000
19 20 21 22 23 24
-18.1875000 -20.1875000 -22.0625000 -20.3125000 -19.1875000 -23.5625000
25 26 27 28 29 30
-25.8333333 -16.9583333 -8.2083333 -3.9583333 -16.4583333 -16.5833333
31 32 33 34 35 36
-14.4583333 -10.4583333 -5.3333333 -7.5833333 -5.4583333 -3.8333333
37 38 39 40 41 42
1.8958333 -7.2291667 -3.4791667 6.7708333 16.2708333 19.1458333
43 44 45 46 47 48
23.2708333 26.2708333 29.3958333 35.1458333 35.2708333 41.8958333
49 50 51 52 53 54
-20.1458333 -26.2708333 -22.5208333 -21.2708333 -19.7708333 -11.8958333
55 56 57 58 59 60
-7.7708333 -3.7708333 0.3541667 -3.8958333 -2.7708333 3.8541667
61 62 63 64 65 66
-0.4166667 -0.5416667 7.2083333 19.4583333 21.9583333 22.8333333
67 68 69 70 71 72
15.9583333 18.9583333 24.0833333 19.8333333 21.9583333 24.5833333
73 74 75 76 77 78
21.3125000 12.1875000 16.9375000 15.1875000 17.6875000 19.5625000
79 80 81 82 83 84
12.6875000 14.6875000 14.8125000 15.5625000 22.6875000 23.3125000
85 86 87 88 89 90
15.0416667 3.9166667 3.6666667 -14.0833333 -16.5833333 -25.7083333
91 92 93 94 95 96
-24.5833333 -31.5833333 -38.4583333 -37.7083333 -44.5833333 -55.9583333
> postscript(file="/var/www/html/rcomp/tmp/6ocae1229456596.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 = 96
Frequency = 1
lag(myerror, k = 1) myerror
0 24.7083333 NA
1 35.5833333 24.7083333
2 21.3333333 35.5833333
3 14.5833333 21.3333333
4 15.0833333 14.5833333
5 14.9583333 15.0833333
6 13.0833333 14.9583333
7 6.0833333 13.0833333
8 -2.7916667 6.0833333
9 -1.0416667 -2.7916667
10 -7.9166667 -1.0416667
11 -10.2916667 -7.9166667
12 -16.5625000 -10.2916667
13 -0.6875000 -16.5625000
14 -14.9375000 -0.6875000
15 -16.6875000 -14.9375000
16 -18.1875000 -16.6875000
17 -22.3125000 -18.1875000
18 -18.1875000 -22.3125000
19 -20.1875000 -18.1875000
20 -22.0625000 -20.1875000
21 -20.3125000 -22.0625000
22 -19.1875000 -20.3125000
23 -23.5625000 -19.1875000
24 -25.8333333 -23.5625000
25 -16.9583333 -25.8333333
26 -8.2083333 -16.9583333
27 -3.9583333 -8.2083333
28 -16.4583333 -3.9583333
29 -16.5833333 -16.4583333
30 -14.4583333 -16.5833333
31 -10.4583333 -14.4583333
32 -5.3333333 -10.4583333
33 -7.5833333 -5.3333333
34 -5.4583333 -7.5833333
35 -3.8333333 -5.4583333
36 1.8958333 -3.8333333
37 -7.2291667 1.8958333
38 -3.4791667 -7.2291667
39 6.7708333 -3.4791667
40 16.2708333 6.7708333
41 19.1458333 16.2708333
42 23.2708333 19.1458333
43 26.2708333 23.2708333
44 29.3958333 26.2708333
45 35.1458333 29.3958333
46 35.2708333 35.1458333
47 41.8958333 35.2708333
48 -20.1458333 41.8958333
49 -26.2708333 -20.1458333
50 -22.5208333 -26.2708333
51 -21.2708333 -22.5208333
52 -19.7708333 -21.2708333
53 -11.8958333 -19.7708333
54 -7.7708333 -11.8958333
55 -3.7708333 -7.7708333
56 0.3541667 -3.7708333
57 -3.8958333 0.3541667
58 -2.7708333 -3.8958333
59 3.8541667 -2.7708333
60 -0.4166667 3.8541667
61 -0.5416667 -0.4166667
62 7.2083333 -0.5416667
63 19.4583333 7.2083333
64 21.9583333 19.4583333
65 22.8333333 21.9583333
66 15.9583333 22.8333333
67 18.9583333 15.9583333
68 24.0833333 18.9583333
69 19.8333333 24.0833333
70 21.9583333 19.8333333
71 24.5833333 21.9583333
72 21.3125000 24.5833333
73 12.1875000 21.3125000
74 16.9375000 12.1875000
75 15.1875000 16.9375000
76 17.6875000 15.1875000
77 19.5625000 17.6875000
78 12.6875000 19.5625000
79 14.6875000 12.6875000
80 14.8125000 14.6875000
81 15.5625000 14.8125000
82 22.6875000 15.5625000
83 23.3125000 22.6875000
84 15.0416667 23.3125000
85 3.9166667 15.0416667
86 3.6666667 3.9166667
87 -14.0833333 3.6666667
88 -16.5833333 -14.0833333
89 -25.7083333 -16.5833333
90 -24.5833333 -25.7083333
91 -31.5833333 -24.5833333
92 -38.4583333 -31.5833333
93 -37.7083333 -38.4583333
94 -44.5833333 -37.7083333
95 -55.9583333 -44.5833333
96 NA -55.9583333
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 35.5833333 24.7083333
[2,] 21.3333333 35.5833333
[3,] 14.5833333 21.3333333
[4,] 15.0833333 14.5833333
[5,] 14.9583333 15.0833333
[6,] 13.0833333 14.9583333
[7,] 6.0833333 13.0833333
[8,] -2.7916667 6.0833333
[9,] -1.0416667 -2.7916667
[10,] -7.9166667 -1.0416667
[11,] -10.2916667 -7.9166667
[12,] -16.5625000 -10.2916667
[13,] -0.6875000 -16.5625000
[14,] -14.9375000 -0.6875000
[15,] -16.6875000 -14.9375000
[16,] -18.1875000 -16.6875000
[17,] -22.3125000 -18.1875000
[18,] -18.1875000 -22.3125000
[19,] -20.1875000 -18.1875000
[20,] -22.0625000 -20.1875000
[21,] -20.3125000 -22.0625000
[22,] -19.1875000 -20.3125000
[23,] -23.5625000 -19.1875000
[24,] -25.8333333 -23.5625000
[25,] -16.9583333 -25.8333333
[26,] -8.2083333 -16.9583333
[27,] -3.9583333 -8.2083333
[28,] -16.4583333 -3.9583333
[29,] -16.5833333 -16.4583333
[30,] -14.4583333 -16.5833333
[31,] -10.4583333 -14.4583333
[32,] -5.3333333 -10.4583333
[33,] -7.5833333 -5.3333333
[34,] -5.4583333 -7.5833333
[35,] -3.8333333 -5.4583333
[36,] 1.8958333 -3.8333333
[37,] -7.2291667 1.8958333
[38,] -3.4791667 -7.2291667
[39,] 6.7708333 -3.4791667
[40,] 16.2708333 6.7708333
[41,] 19.1458333 16.2708333
[42,] 23.2708333 19.1458333
[43,] 26.2708333 23.2708333
[44,] 29.3958333 26.2708333
[45,] 35.1458333 29.3958333
[46,] 35.2708333 35.1458333
[47,] 41.8958333 35.2708333
[48,] -20.1458333 41.8958333
[49,] -26.2708333 -20.1458333
[50,] -22.5208333 -26.2708333
[51,] -21.2708333 -22.5208333
[52,] -19.7708333 -21.2708333
[53,] -11.8958333 -19.7708333
[54,] -7.7708333 -11.8958333
[55,] -3.7708333 -7.7708333
[56,] 0.3541667 -3.7708333
[57,] -3.8958333 0.3541667
[58,] -2.7708333 -3.8958333
[59,] 3.8541667 -2.7708333
[60,] -0.4166667 3.8541667
[61,] -0.5416667 -0.4166667
[62,] 7.2083333 -0.5416667
[63,] 19.4583333 7.2083333
[64,] 21.9583333 19.4583333
[65,] 22.8333333 21.9583333
[66,] 15.9583333 22.8333333
[67,] 18.9583333 15.9583333
[68,] 24.0833333 18.9583333
[69,] 19.8333333 24.0833333
[70,] 21.9583333 19.8333333
[71,] 24.5833333 21.9583333
[72,] 21.3125000 24.5833333
[73,] 12.1875000 21.3125000
[74,] 16.9375000 12.1875000
[75,] 15.1875000 16.9375000
[76,] 17.6875000 15.1875000
[77,] 19.5625000 17.6875000
[78,] 12.6875000 19.5625000
[79,] 14.6875000 12.6875000
[80,] 14.8125000 14.6875000
[81,] 15.5625000 14.8125000
[82,] 22.6875000 15.5625000
[83,] 23.3125000 22.6875000
[84,] 15.0416667 23.3125000
[85,] 3.9166667 15.0416667
[86,] 3.6666667 3.9166667
[87,] -14.0833333 3.6666667
[88,] -16.5833333 -14.0833333
[89,] -25.7083333 -16.5833333
[90,] -24.5833333 -25.7083333
[91,] -31.5833333 -24.5833333
[92,] -38.4583333 -31.5833333
[93,] -37.7083333 -38.4583333
[94,] -44.5833333 -37.7083333
[95,] -55.9583333 -44.5833333
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 35.5833333 24.7083333
2 21.3333333 35.5833333
3 14.5833333 21.3333333
4 15.0833333 14.5833333
5 14.9583333 15.0833333
6 13.0833333 14.9583333
7 6.0833333 13.0833333
8 -2.7916667 6.0833333
9 -1.0416667 -2.7916667
10 -7.9166667 -1.0416667
11 -10.2916667 -7.9166667
12 -16.5625000 -10.2916667
13 -0.6875000 -16.5625000
14 -14.9375000 -0.6875000
15 -16.6875000 -14.9375000
16 -18.1875000 -16.6875000
17 -22.3125000 -18.1875000
18 -18.1875000 -22.3125000
19 -20.1875000 -18.1875000
20 -22.0625000 -20.1875000
21 -20.3125000 -22.0625000
22 -19.1875000 -20.3125000
23 -23.5625000 -19.1875000
24 -25.8333333 -23.5625000
25 -16.9583333 -25.8333333
26 -8.2083333 -16.9583333
27 -3.9583333 -8.2083333
28 -16.4583333 -3.9583333
29 -16.5833333 -16.4583333
30 -14.4583333 -16.5833333
31 -10.4583333 -14.4583333
32 -5.3333333 -10.4583333
33 -7.5833333 -5.3333333
34 -5.4583333 -7.5833333
35 -3.8333333 -5.4583333
36 1.8958333 -3.8333333
37 -7.2291667 1.8958333
38 -3.4791667 -7.2291667
39 6.7708333 -3.4791667
40 16.2708333 6.7708333
41 19.1458333 16.2708333
42 23.2708333 19.1458333
43 26.2708333 23.2708333
44 29.3958333 26.2708333
45 35.1458333 29.3958333
46 35.2708333 35.1458333
47 41.8958333 35.2708333
48 -20.1458333 41.8958333
49 -26.2708333 -20.1458333
50 -22.5208333 -26.2708333
51 -21.2708333 -22.5208333
52 -19.7708333 -21.2708333
53 -11.8958333 -19.7708333
54 -7.7708333 -11.8958333
55 -3.7708333 -7.7708333
56 0.3541667 -3.7708333
57 -3.8958333 0.3541667
58 -2.7708333 -3.8958333
59 3.8541667 -2.7708333
60 -0.4166667 3.8541667
61 -0.5416667 -0.4166667
62 7.2083333 -0.5416667
63 19.4583333 7.2083333
64 21.9583333 19.4583333
65 22.8333333 21.9583333
66 15.9583333 22.8333333
67 18.9583333 15.9583333
68 24.0833333 18.9583333
69 19.8333333 24.0833333
70 21.9583333 19.8333333
71 24.5833333 21.9583333
72 21.3125000 24.5833333
73 12.1875000 21.3125000
74 16.9375000 12.1875000
75 15.1875000 16.9375000
76 17.6875000 15.1875000
77 19.5625000 17.6875000
78 12.6875000 19.5625000
79 14.6875000 12.6875000
80 14.8125000 14.6875000
81 15.5625000 14.8125000
82 22.6875000 15.5625000
83 23.3125000 22.6875000
84 15.0416667 23.3125000
85 3.9166667 15.0416667
86 3.6666667 3.9166667
87 -14.0833333 3.6666667
88 -16.5833333 -14.0833333
89 -25.7083333 -16.5833333
90 -24.5833333 -25.7083333
91 -31.5833333 -24.5833333
92 -38.4583333 -31.5833333
93 -37.7083333 -38.4583333
94 -44.5833333 -37.7083333
95 -55.9583333 -44.5833333
> 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/75qnt1229456596.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/8qyto1229456596.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/95d821229456596.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/10tke31229456596.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/11ax8e1229456596.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/12fjy11229456596.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/133dxx1229456596.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/14kcen1229456596.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/153viu1229456596.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/169xkb1229456596.tab")
+ }
>
> system("convert tmp/15muw1229456596.ps tmp/15muw1229456596.png")
> system("convert tmp/26jle1229456596.ps tmp/26jle1229456596.png")
> system("convert tmp/399np1229456596.ps tmp/399np1229456596.png")
> system("convert tmp/45kyj1229456596.ps tmp/45kyj1229456596.png")
> system("convert tmp/5uaxl1229456596.ps tmp/5uaxl1229456596.png")
> system("convert tmp/6ocae1229456596.ps tmp/6ocae1229456596.png")
> system("convert tmp/75qnt1229456596.ps tmp/75qnt1229456596.png")
> system("convert tmp/8qyto1229456596.ps tmp/8qyto1229456596.png")
> system("convert tmp/95d821229456596.ps tmp/95d821229456596.png")
> system("convert tmp/10tke31229456596.ps tmp/10tke31229456596.png")
>
>
> proc.time()
user system elapsed
2.977 1.654 3.439