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,1,481,1,462,1,457,1,442,1,439,1,488,1,521,1,501,1,485,1,464,1,460,1,467,1,460,1,448,1,443,1,436,1,431,1,484,1,510,1,513,1,503,1,471,1,471,1,476,1,475,1,470,1,461,1,455,1,456,1,517,1,525,1,523,1,519,1,509,1,512,1,519,0,517,0,510,0,509,0,501,0,507,0,569,0,580,0,578,0,565,0,547,0,555,0,562,0,561,0,555,0,544,0,537,0,543,0,594,0,611,0,613,0,611,0,594,0,595,0,591,0,589,0,584,0,573,0,567,0,569,0,621,0,629,0,628,0,612,0,595,0,597,0,593,0,590,0,580,0,574,0,573,0,573,0,620,0,626,0,620,0,588,0,566,0,557,0,561,0,549,0,532,0,526,0,511,0,499,0,555,0,565,0,542,0,527,0,510,0,514,0,517,0,508,0,493,0,490,0,469,0,478,0),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 = 'No Linear Trend'
> par2 = 'Do not include Seasonal Dummies'
> par1 = '1'
> #'GNU S' R Code compiled by R2WASP v. 1.0.44 ()
> #Author: Prof. Dr. P. Wessa
> #To cite this work: AUTHOR(S), (YEAR), YOUR SOFTWARE TITLE (vNUMBER) in Free Statistics Software (v$_version), Office for Research Development and Education, URL http://www.wessa.net/rwasp_YOURPAGE.wasp/
> #Source of accompanying publication: Office for Research, Development, and Education
> #Technical description: Write here your technical program description (don't use hard returns!)
> library(lattice)
> library(lmtest)
Loading required package: zoo
Attaching package: 'zoo'
The following object(s) are masked from package:base :
as.Date.numeric
> n25 <- 25 #minimum number of obs. for Goldfeld-Quandt test
> par1 <- as.numeric(par1)
> x <- t(y)
> k <- length(x[1,])
> n <- length(x[,1])
> x1 <- cbind(x[,par1], x[,1:k!=par1])
> mycolnames <- c(colnames(x)[par1], colnames(x)[1:k!=par1])
> colnames(x1) <- mycolnames #colnames(x)[par1]
> x <- x1
> if (par3 == 'First Differences'){
+ x2 <- array(0, dim=c(n-1,k), dimnames=list(1:(n-1), paste('(1-B)',colnames(x),sep='')))
+ for (i in 1:n-1) {
+ for (j in 1:k) {
+ x2[i,j] <- x[i+1,j] - x[i,j]
+ }
+ }
+ x <- x2
+ }
> if (par2 == 'Include Monthly Dummies'){
+ x2 <- array(0, dim=c(n,11), dimnames=list(1:n, paste('M', seq(1:11), sep ='')))
+ for (i in 1:11){
+ x2[seq(i,n,12),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> if (par2 == 'Include Quarterly Dummies'){
+ x2 <- array(0, dim=c(n,3), dimnames=list(1:n, paste('Q', seq(1:3), sep ='')))
+ for (i in 1:3){
+ x2[seq(i,n,4),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> k <- length(x[1,])
> if (par3 == 'Linear Trend'){
+ x <- cbind(x, c(1:n))
+ colnames(x)[k+1] <- 't'
+ }
> x
Aantal_werklozen_(*1000) dummyvariabele
1 493 1
2 481 1
3 462 1
4 457 1
5 442 1
6 439 1
7 488 1
8 521 1
9 501 1
10 485 1
11 464 1
12 460 1
13 467 1
14 460 1
15 448 1
16 443 1
17 436 1
18 431 1
19 484 1
20 510 1
21 513 1
22 503 1
23 471 1
24 471 1
25 476 1
26 475 1
27 470 1
28 461 1
29 455 1
30 456 1
31 517 1
32 525 1
33 523 1
34 519 1
35 509 1
36 512 1
37 519 0
38 517 0
39 510 0
40 509 0
41 501 0
42 507 0
43 569 0
44 580 0
45 578 0
46 565 0
47 547 0
48 555 0
49 562 0
50 561 0
51 555 0
52 544 0
53 537 0
54 543 0
55 594 0
56 611 0
57 613 0
58 611 0
59 594 0
60 595 0
61 591 0
62 589 0
63 584 0
64 573 0
65 567 0
66 569 0
67 621 0
68 629 0
69 628 0
70 612 0
71 595 0
72 597 0
73 593 0
74 590 0
75 580 0
76 574 0
77 573 0
78 573 0
79 620 0
80 626 0
81 620 0
82 588 0
83 566 0
84 557 0
85 561 0
86 549 0
87 532 0
88 526 0
89 511 0
90 499 0
91 555 0
92 565 0
93 542 0
94 527 0
95 510 0
96 514 0
97 517 0
98 508 0
99 493 0
100 490 0
101 469 0
102 478 0
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) dummyvariabele
559.67 -81.11
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-90.667 -23.333 1.833 30.417 69.333
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 559.667 4.530 123.54 <2e-16 ***
dummyvariabele -81.111 7.626 -10.64 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 36.8 on 100 degrees of freedom
Multiple R-squared: 0.5308, Adjusted R-squared: 0.5261
F-statistic: 113.1 on 1 and 100 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.2366864490 0.4733728981 0.76331355
[2,] 0.1970142797 0.3940285594 0.80298572
[3,] 0.1540188655 0.3080377309 0.84598113
[4,] 0.3092223414 0.6184446827 0.69077766
[5,] 0.2642438687 0.5284877374 0.73575613
[6,] 0.1808385113 0.3616770226 0.81916149
[7,] 0.1236226426 0.2472452852 0.87637736
[8,] 0.0856643050 0.1713286100 0.91433570
[9,] 0.0529874422 0.1059748844 0.94701256
[10,] 0.0344408980 0.0688817960 0.96555910
[11,] 0.0286414084 0.0572828167 0.97135859
[12,] 0.0267438157 0.0534876314 0.97325618
[13,] 0.0302277221 0.0604554442 0.96977228
[14,] 0.0382548125 0.0765096250 0.96174519
[15,] 0.0276044748 0.0552089496 0.97239553
[16,] 0.0379920563 0.0759841126 0.96200794
[17,] 0.0503223251 0.1006446502 0.94967767
[18,] 0.0471703743 0.0943407487 0.95282963
[19,] 0.0316535074 0.0633070148 0.96834649
[20,] 0.0208133811 0.0416267622 0.97918662
[21,] 0.0133672136 0.0267344272 0.98663279
[22,] 0.0084159773 0.0168319546 0.99158402
[23,] 0.0052899097 0.0105798194 0.99471009
[24,] 0.0036546298 0.0073092597 0.99634537
[25,] 0.0029197267 0.0058394534 0.99708027
[26,] 0.0024234205 0.0048468410 0.99757658
[27,] 0.0038425096 0.0076850191 0.99615749
[28,] 0.0071247458 0.0142494917 0.99287525
[29,] 0.0103530081 0.0207060163 0.98964699
[30,] 0.0121233932 0.0242467865 0.98787661
[31,] 0.0107988651 0.0215977303 0.98920113
[32,] 0.0099427892 0.0198855784 0.99005721
[33,] 0.0073495610 0.0146991220 0.99265044
[34,] 0.0055098040 0.0110196080 0.99449020
[35,] 0.0044543152 0.0089086304 0.99554568
[36,] 0.0036730332 0.0073460664 0.99632697
[37,] 0.0034987000 0.0069974001 0.99650130
[38,] 0.0030781884 0.0061563768 0.99692181
[39,] 0.0050900125 0.0101800249 0.99490999
[40,] 0.0085834092 0.0171668183 0.99141659
[41,] 0.0105796532 0.0211593063 0.98942035
[42,] 0.0089777357 0.0179554713 0.99102226
[43,] 0.0063324119 0.0126648237 0.99366759
[44,] 0.0045322683 0.0090645366 0.99546773
[45,] 0.0033640636 0.0067281272 0.99663594
[46,] 0.0023988656 0.0047977312 0.99760113
[47,] 0.0015972786 0.0031945571 0.99840272
[48,] 0.0010416002 0.0020832004 0.99895840
[49,] 0.0007075039 0.0014150078 0.99929250
[50,] 0.0004542037 0.0009084075 0.99954580
[51,] 0.0006450646 0.0012901293 0.99935494
[52,] 0.0015620743 0.0031241486 0.99843793
[53,] 0.0033214406 0.0066428812 0.99667856
[54,] 0.0056516010 0.0113032019 0.99434840
[55,] 0.0056110685 0.0112221371 0.99438893
[56,] 0.0055640400 0.0111280800 0.99443596
[57,] 0.0049922036 0.0099844073 0.99500780
[58,] 0.0042605438 0.0085210876 0.99573946
[59,] 0.0033107485 0.0066214971 0.99668925
[60,] 0.0022309529 0.0044619059 0.99776905
[61,] 0.0014197334 0.0028394668 0.99858027
[62,] 0.0008958826 0.0017917651 0.99910412
[63,] 0.0019577097 0.0039154193 0.99804229
[64,] 0.0055477830 0.0110955660 0.99445222
[65,] 0.0139879760 0.0279759519 0.98601202
[66,] 0.0205536566 0.0411073131 0.97944634
[67,] 0.0207714012 0.0415428025 0.97922860
[68,] 0.0223696899 0.0447393799 0.97763031
[69,] 0.0228753753 0.0457507505 0.97712462
[70,] 0.0227391679 0.0454783359 0.97726083
[71,] 0.0196468189 0.0392936379 0.98035318
[72,] 0.0158821575 0.0317643149 0.98411784
[73,] 0.0127748203 0.0255496407 0.98722518
[74,] 0.0103725707 0.0207451414 0.98962743
[75,] 0.0329632513 0.0659265025 0.96703675
[76,] 0.1372425126 0.2744850252 0.86275749
[77,] 0.4162015112 0.8324030223 0.58379849
[78,] 0.5753539319 0.8492921362 0.42464607
[79,] 0.6246268602 0.7507462796 0.37537314
[80,] 0.6452464934 0.7095070132 0.35475351
[81,] 0.7035462711 0.5929074577 0.29645373
[82,] 0.7175328556 0.5649342887 0.28246714
[83,] 0.6820021399 0.6359957201 0.31799786
[84,] 0.6331719702 0.7336560595 0.36682803
[85,] 0.5782824080 0.8434351840 0.42171759
[86,] 0.5456082656 0.9087834688 0.45439173
[87,] 0.6099099673 0.7801800655 0.39009003
[88,] 0.8303644312 0.3392711376 0.16963557
[89,] 0.8994031644 0.2011936711 0.10059684
[90,] 0.9125520271 0.1748959458 0.08744797
[91,] 0.8690994992 0.2618010016 0.13090050
[92,] 0.8340777153 0.3318445695 0.16592228
[93,] 0.8508416924 0.2983166152 0.14915831
> postscript(file="/var/www/html/rcomp/tmp/1mkry1228657795.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/22enz1228657795.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/3cnly1228657795.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/4cn8e1228657795.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/5vs311228657795.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 7
14.444444 2.444444 -16.555556 -21.555556 -36.555556 -39.555556 9.444444
8 9 10 11 12 13 14
42.444444 22.444444 6.444444 -14.555556 -18.555556 -11.555556 -18.555556
15 16 17 18 19 20 21
-30.555556 -35.555556 -42.555556 -47.555556 5.444444 31.444444 34.444444
22 23 24 25 26 27 28
24.444444 -7.555556 -7.555556 -2.555556 -3.555556 -8.555556 -17.555556
29 30 31 32 33 34 35
-23.555556 -22.555556 38.444444 46.444444 44.444444 40.444444 30.444444
36 37 38 39 40 41 42
33.444444 -40.666667 -42.666667 -49.666667 -50.666667 -58.666667 -52.666667
43 44 45 46 47 48 49
9.333333 20.333333 18.333333 5.333333 -12.666667 -4.666667 2.333333
50 51 52 53 54 55 56
1.333333 -4.666667 -15.666667 -22.666667 -16.666667 34.333333 51.333333
57 58 59 60 61 62 63
53.333333 51.333333 34.333333 35.333333 31.333333 29.333333 24.333333
64 65 66 67 68 69 70
13.333333 7.333333 9.333333 61.333333 69.333333 68.333333 52.333333
71 72 73 74 75 76 77
35.333333 37.333333 33.333333 30.333333 20.333333 14.333333 13.333333
78 79 80 81 82 83 84
13.333333 60.333333 66.333333 60.333333 28.333333 6.333333 -2.666667
85 86 87 88 89 90 91
1.333333 -10.666667 -27.666667 -33.666667 -48.666667 -60.666667 -4.666667
92 93 94 95 96 97 98
5.333333 -17.666667 -32.666667 -49.666667 -45.666667 -42.666667 -51.666667
99 100 101 102
-66.666667 -69.666667 -90.666667 -81.666667
> postscript(file="/var/www/html/rcomp/tmp/648f51228657795.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 14.444444 NA
1 2.444444 14.444444
2 -16.555556 2.444444
3 -21.555556 -16.555556
4 -36.555556 -21.555556
5 -39.555556 -36.555556
6 9.444444 -39.555556
7 42.444444 9.444444
8 22.444444 42.444444
9 6.444444 22.444444
10 -14.555556 6.444444
11 -18.555556 -14.555556
12 -11.555556 -18.555556
13 -18.555556 -11.555556
14 -30.555556 -18.555556
15 -35.555556 -30.555556
16 -42.555556 -35.555556
17 -47.555556 -42.555556
18 5.444444 -47.555556
19 31.444444 5.444444
20 34.444444 31.444444
21 24.444444 34.444444
22 -7.555556 24.444444
23 -7.555556 -7.555556
24 -2.555556 -7.555556
25 -3.555556 -2.555556
26 -8.555556 -3.555556
27 -17.555556 -8.555556
28 -23.555556 -17.555556
29 -22.555556 -23.555556
30 38.444444 -22.555556
31 46.444444 38.444444
32 44.444444 46.444444
33 40.444444 44.444444
34 30.444444 40.444444
35 33.444444 30.444444
36 -40.666667 33.444444
37 -42.666667 -40.666667
38 -49.666667 -42.666667
39 -50.666667 -49.666667
40 -58.666667 -50.666667
41 -52.666667 -58.666667
42 9.333333 -52.666667
43 20.333333 9.333333
44 18.333333 20.333333
45 5.333333 18.333333
46 -12.666667 5.333333
47 -4.666667 -12.666667
48 2.333333 -4.666667
49 1.333333 2.333333
50 -4.666667 1.333333
51 -15.666667 -4.666667
52 -22.666667 -15.666667
53 -16.666667 -22.666667
54 34.333333 -16.666667
55 51.333333 34.333333
56 53.333333 51.333333
57 51.333333 53.333333
58 34.333333 51.333333
59 35.333333 34.333333
60 31.333333 35.333333
61 29.333333 31.333333
62 24.333333 29.333333
63 13.333333 24.333333
64 7.333333 13.333333
65 9.333333 7.333333
66 61.333333 9.333333
67 69.333333 61.333333
68 68.333333 69.333333
69 52.333333 68.333333
70 35.333333 52.333333
71 37.333333 35.333333
72 33.333333 37.333333
73 30.333333 33.333333
74 20.333333 30.333333
75 14.333333 20.333333
76 13.333333 14.333333
77 13.333333 13.333333
78 60.333333 13.333333
79 66.333333 60.333333
80 60.333333 66.333333
81 28.333333 60.333333
82 6.333333 28.333333
83 -2.666667 6.333333
84 1.333333 -2.666667
85 -10.666667 1.333333
86 -27.666667 -10.666667
87 -33.666667 -27.666667
88 -48.666667 -33.666667
89 -60.666667 -48.666667
90 -4.666667 -60.666667
91 5.333333 -4.666667
92 -17.666667 5.333333
93 -32.666667 -17.666667
94 -49.666667 -32.666667
95 -45.666667 -49.666667
96 -42.666667 -45.666667
97 -51.666667 -42.666667
98 -66.666667 -51.666667
99 -69.666667 -66.666667
100 -90.666667 -69.666667
101 -81.666667 -90.666667
102 NA -81.666667
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 2.444444 14.444444
[2,] -16.555556 2.444444
[3,] -21.555556 -16.555556
[4,] -36.555556 -21.555556
[5,] -39.555556 -36.555556
[6,] 9.444444 -39.555556
[7,] 42.444444 9.444444
[8,] 22.444444 42.444444
[9,] 6.444444 22.444444
[10,] -14.555556 6.444444
[11,] -18.555556 -14.555556
[12,] -11.555556 -18.555556
[13,] -18.555556 -11.555556
[14,] -30.555556 -18.555556
[15,] -35.555556 -30.555556
[16,] -42.555556 -35.555556
[17,] -47.555556 -42.555556
[18,] 5.444444 -47.555556
[19,] 31.444444 5.444444
[20,] 34.444444 31.444444
[21,] 24.444444 34.444444
[22,] -7.555556 24.444444
[23,] -7.555556 -7.555556
[24,] -2.555556 -7.555556
[25,] -3.555556 -2.555556
[26,] -8.555556 -3.555556
[27,] -17.555556 -8.555556
[28,] -23.555556 -17.555556
[29,] -22.555556 -23.555556
[30,] 38.444444 -22.555556
[31,] 46.444444 38.444444
[32,] 44.444444 46.444444
[33,] 40.444444 44.444444
[34,] 30.444444 40.444444
[35,] 33.444444 30.444444
[36,] -40.666667 33.444444
[37,] -42.666667 -40.666667
[38,] -49.666667 -42.666667
[39,] -50.666667 -49.666667
[40,] -58.666667 -50.666667
[41,] -52.666667 -58.666667
[42,] 9.333333 -52.666667
[43,] 20.333333 9.333333
[44,] 18.333333 20.333333
[45,] 5.333333 18.333333
[46,] -12.666667 5.333333
[47,] -4.666667 -12.666667
[48,] 2.333333 -4.666667
[49,] 1.333333 2.333333
[50,] -4.666667 1.333333
[51,] -15.666667 -4.666667
[52,] -22.666667 -15.666667
[53,] -16.666667 -22.666667
[54,] 34.333333 -16.666667
[55,] 51.333333 34.333333
[56,] 53.333333 51.333333
[57,] 51.333333 53.333333
[58,] 34.333333 51.333333
[59,] 35.333333 34.333333
[60,] 31.333333 35.333333
[61,] 29.333333 31.333333
[62,] 24.333333 29.333333
[63,] 13.333333 24.333333
[64,] 7.333333 13.333333
[65,] 9.333333 7.333333
[66,] 61.333333 9.333333
[67,] 69.333333 61.333333
[68,] 68.333333 69.333333
[69,] 52.333333 68.333333
[70,] 35.333333 52.333333
[71,] 37.333333 35.333333
[72,] 33.333333 37.333333
[73,] 30.333333 33.333333
[74,] 20.333333 30.333333
[75,] 14.333333 20.333333
[76,] 13.333333 14.333333
[77,] 13.333333 13.333333
[78,] 60.333333 13.333333
[79,] 66.333333 60.333333
[80,] 60.333333 66.333333
[81,] 28.333333 60.333333
[82,] 6.333333 28.333333
[83,] -2.666667 6.333333
[84,] 1.333333 -2.666667
[85,] -10.666667 1.333333
[86,] -27.666667 -10.666667
[87,] -33.666667 -27.666667
[88,] -48.666667 -33.666667
[89,] -60.666667 -48.666667
[90,] -4.666667 -60.666667
[91,] 5.333333 -4.666667
[92,] -17.666667 5.333333
[93,] -32.666667 -17.666667
[94,] -49.666667 -32.666667
[95,] -45.666667 -49.666667
[96,] -42.666667 -45.666667
[97,] -51.666667 -42.666667
[98,] -66.666667 -51.666667
[99,] -69.666667 -66.666667
[100,] -90.666667 -69.666667
[101,] -81.666667 -90.666667
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 2.444444 14.444444
2 -16.555556 2.444444
3 -21.555556 -16.555556
4 -36.555556 -21.555556
5 -39.555556 -36.555556
6 9.444444 -39.555556
7 42.444444 9.444444
8 22.444444 42.444444
9 6.444444 22.444444
10 -14.555556 6.444444
11 -18.555556 -14.555556
12 -11.555556 -18.555556
13 -18.555556 -11.555556
14 -30.555556 -18.555556
15 -35.555556 -30.555556
16 -42.555556 -35.555556
17 -47.555556 -42.555556
18 5.444444 -47.555556
19 31.444444 5.444444
20 34.444444 31.444444
21 24.444444 34.444444
22 -7.555556 24.444444
23 -7.555556 -7.555556
24 -2.555556 -7.555556
25 -3.555556 -2.555556
26 -8.555556 -3.555556
27 -17.555556 -8.555556
28 -23.555556 -17.555556
29 -22.555556 -23.555556
30 38.444444 -22.555556
31 46.444444 38.444444
32 44.444444 46.444444
33 40.444444 44.444444
34 30.444444 40.444444
35 33.444444 30.444444
36 -40.666667 33.444444
37 -42.666667 -40.666667
38 -49.666667 -42.666667
39 -50.666667 -49.666667
40 -58.666667 -50.666667
41 -52.666667 -58.666667
42 9.333333 -52.666667
43 20.333333 9.333333
44 18.333333 20.333333
45 5.333333 18.333333
46 -12.666667 5.333333
47 -4.666667 -12.666667
48 2.333333 -4.666667
49 1.333333 2.333333
50 -4.666667 1.333333
51 -15.666667 -4.666667
52 -22.666667 -15.666667
53 -16.666667 -22.666667
54 34.333333 -16.666667
55 51.333333 34.333333
56 53.333333 51.333333
57 51.333333 53.333333
58 34.333333 51.333333
59 35.333333 34.333333
60 31.333333 35.333333
61 29.333333 31.333333
62 24.333333 29.333333
63 13.333333 24.333333
64 7.333333 13.333333
65 9.333333 7.333333
66 61.333333 9.333333
67 69.333333 61.333333
68 68.333333 69.333333
69 52.333333 68.333333
70 35.333333 52.333333
71 37.333333 35.333333
72 33.333333 37.333333
73 30.333333 33.333333
74 20.333333 30.333333
75 14.333333 20.333333
76 13.333333 14.333333
77 13.333333 13.333333
78 60.333333 13.333333
79 66.333333 60.333333
80 60.333333 66.333333
81 28.333333 60.333333
82 6.333333 28.333333
83 -2.666667 6.333333
84 1.333333 -2.666667
85 -10.666667 1.333333
86 -27.666667 -10.666667
87 -33.666667 -27.666667
88 -48.666667 -33.666667
89 -60.666667 -48.666667
90 -4.666667 -60.666667
91 5.333333 -4.666667
92 -17.666667 5.333333
93 -32.666667 -17.666667
94 -49.666667 -32.666667
95 -45.666667 -49.666667
96 -42.666667 -45.666667
97 -51.666667 -42.666667
98 -66.666667 -51.666667
99 -69.666667 -66.666667
100 -90.666667 -69.666667
101 -81.666667 -90.666667
> 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/7yoox1228657795.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/8oxcx1228657795.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/9gyg61228657796.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/101nci1228657796.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/119wb31228657796.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/1201fy1228657796.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/13z7n01228657796.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/142gxm1228657796.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/15al3h1228657796.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/16m1vb1228657796.tab")
+ }
>
> system("convert tmp/1mkry1228657795.ps tmp/1mkry1228657795.png")
> system("convert tmp/22enz1228657795.ps tmp/22enz1228657795.png")
> system("convert tmp/3cnly1228657795.ps tmp/3cnly1228657795.png")
> system("convert tmp/4cn8e1228657795.ps tmp/4cn8e1228657795.png")
> system("convert tmp/5vs311228657795.ps tmp/5vs311228657795.png")
> system("convert tmp/648f51228657795.ps tmp/648f51228657795.png")
> system("convert tmp/7yoox1228657795.ps tmp/7yoox1228657795.png")
> system("convert tmp/8oxcx1228657795.ps tmp/8oxcx1228657795.png")
> system("convert tmp/9gyg61228657796.ps tmp/9gyg61228657796.png")
> system("convert tmp/101nci1228657796.ps tmp/101nci1228657796.png")
>
>
> proc.time()
user system elapsed
2.970 1.637 3.875