R version 2.13.0 (2011-04-13)
Copyright (C) 2011 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: i486-pc-linux-gnu (32-bit)
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(72772
+ ,26073
+ ,22274
+ ,45104
+ ,18103
+ ,14819
+ ,44525
+ ,15100
+ ,15136
+ ,41169
+ ,14738
+ ,13704
+ ,31118
+ ,22259
+ ,19638
+ ,28435
+ ,10277
+ ,7551
+ ,22162
+ ,6225
+ ,8019
+ ,20202
+ ,7663
+ ,6509
+ ,17773
+ ,6618
+ ,6634
+ ,17094
+ ,9945
+ ,11166
+ ,15153
+ ,7590
+ ,7508
+ ,11218
+ ,4293
+ ,4275
+ ,10796
+ ,4656
+ ,4944
+ ,9594
+ ,5145
+ ,5441
+ ,9309
+ ,2001
+ ,1689
+ ,8556
+ ,1779
+ ,1522
+ ,8041
+ ,1609
+ ,1416
+ ,7639
+ ,2191
+ ,1594
+ ,6884
+ ,1617
+ ,1909
+ ,6642
+ ,2554
+ ,2599
+ ,6321
+ ,2198
+ ,1262
+ ,6216
+ ,1578
+ ,1199
+ ,5865
+ ,3446
+ ,4404
+ ,5799
+ ,1380
+ ,1166
+ ,5695
+ ,1249
+ ,1122
+ ,5644
+ ,1223
+ ,886
+ ,5446
+ ,834
+ ,778
+ ,5395
+ ,3754
+ ,4436
+ ,5363
+ ,2283
+ ,1890
+ ,5338
+ ,3028
+ ,3107
+ ,5160
+ ,1100
+ ,1038
+ ,5091
+ ,457
+ ,300
+ ,5057
+ ,1201
+ ,988
+ ,5039
+ ,2192
+ ,2008
+ ,4880
+ ,1508
+ ,1522
+ ,4735
+ ,1393
+ ,1336
+ ,4693
+ ,952
+ ,976
+ ,4653
+ ,1032
+ ,798
+ ,4586
+ ,1279
+ ,869
+ ,4398
+ ,1370
+ ,1260
+ ,3974
+ ,649
+ ,578
+ ,3858
+ ,1900
+ ,2359
+ ,3826
+ ,666
+ ,736
+ ,3819
+ ,1313
+ ,1690
+ ,3556
+ ,1353
+ ,1201
+ ,3372
+ ,1500
+ ,813
+ ,3193
+ ,877
+ ,778
+ ,3126
+ ,874
+ ,687
+ ,3104
+ ,1133
+ ,1270
+ ,2967
+ ,754
+ ,671
+ ,2848
+ ,695
+ ,1559
+ ,2748
+ ,609
+ ,489
+ ,2649
+ ,696
+ ,773
+ ,2625
+ ,756
+ ,629
+ ,2572
+ ,670
+ ,637
+ ,2548
+ ,301
+ ,277
+ ,2477
+ ,630
+ ,776
+ ,2442
+ ,798
+ ,1651
+ ,2392
+ ,436
+ ,377
+ ,2372
+ ,388
+ ,222
+ ,2346
+ ,864
+ ,1068
+ ,2251
+ ,497
+ ,399
+ ,2230
+ ,449
+ ,547
+ ,2225
+ ,919
+ ,668
+ ,2220
+ ,536
+ ,451
+ ,2205
+ ,673
+ ,724
+ ,2193
+ ,837
+ ,853
+ ,2116
+ ,534
+ ,434
+ ,2102
+ ,845
+ ,730
+ ,2099
+ ,626
+ ,612)
+ ,dim=c(3
+ ,70)
+ ,dimnames=list(c('weekdag'
+ ,'zaterdag'
+ ,'zondag')
+ ,1:70))
> y <- array(NA,dim=c(3,70),dimnames=list(c('weekdag','zaterdag','zondag'),1:70))
> 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 = '3'
> #'GNU S' R Code compiled by R2WASP v. 1.0.44 ()
> #Author: Prof. Dr. P. Wessa
> #To cite this work: AUTHOR(S), (YEAR), YOUR SOFTWARE TITLE (vNUMBER) in Free Statistics Software (v$_version), Office for Research Development and Education, URL http://www.wessa.net/rwasp_YOURPAGE.wasp/
> #Source of accompanying publication: Office for Research, Development, and Education
> #Technical description: Write here your technical program description (don't use hard returns!)
> library(lattice)
> library(lmtest)
Loading required package: zoo
> n25 <- 25 #minimum number of obs. for Goldfeld-Quandt test
> par1 <- as.numeric(par1)
> x <- t(y)
> k <- length(x[1,])
> n <- length(x[,1])
> x1 <- cbind(x[,par1], x[,1:k!=par1])
> mycolnames <- c(colnames(x)[par1], colnames(x)[1:k!=par1])
> colnames(x1) <- mycolnames #colnames(x)[par1]
> x <- x1
> if (par3 == 'First Differences'){
+ x2 <- array(0, dim=c(n-1,k), dimnames=list(1:(n-1), paste('(1-B)',colnames(x),sep='')))
+ for (i in 1:n-1) {
+ for (j in 1:k) {
+ x2[i,j] <- x[i+1,j] - x[i,j]
+ }
+ }
+ x <- x2
+ }
> if (par2 == 'Include Monthly Dummies'){
+ x2 <- array(0, dim=c(n,11), dimnames=list(1:n, paste('M', seq(1:11), sep ='')))
+ for (i in 1:11){
+ x2[seq(i,n,12),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> if (par2 == 'Include Quarterly Dummies'){
+ x2 <- array(0, dim=c(n,3), dimnames=list(1:n, paste('Q', seq(1:3), sep ='')))
+ for (i in 1:3){
+ x2[seq(i,n,4),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> k <- length(x[1,])
> if (par3 == 'Linear Trend'){
+ x <- cbind(x, c(1:n))
+ colnames(x)[k+1] <- 't'
+ }
> x
zondag weekdag zaterdag
1 22274 72772 26073
2 14819 45104 18103
3 15136 44525 15100
4 13704 41169 14738
5 19638 31118 22259
6 7551 28435 10277
7 8019 22162 6225
8 6509 20202 7663
9 6634 17773 6618
10 11166 17094 9945
11 7508 15153 7590
12 4275 11218 4293
13 4944 10796 4656
14 5441 9594 5145
15 1689 9309 2001
16 1522 8556 1779
17 1416 8041 1609
18 1594 7639 2191
19 1909 6884 1617
20 2599 6642 2554
21 1262 6321 2198
22 1199 6216 1578
23 4404 5865 3446
24 1166 5799 1380
25 1122 5695 1249
26 886 5644 1223
27 778 5446 834
28 4436 5395 3754
29 1890 5363 2283
30 3107 5338 3028
31 1038 5160 1100
32 300 5091 457
33 988 5057 1201
34 2008 5039 2192
35 1522 4880 1508
36 1336 4735 1393
37 976 4693 952
38 798 4653 1032
39 869 4586 1279
40 1260 4398 1370
41 578 3974 649
42 2359 3858 1900
43 736 3826 666
44 1690 3819 1313
45 1201 3556 1353
46 813 3372 1500
47 778 3193 877
48 687 3126 874
49 1270 3104 1133
50 671 2967 754
51 1559 2848 695
52 489 2748 609
53 773 2649 696
54 629 2625 756
55 637 2572 670
56 277 2548 301
57 776 2477 630
58 1651 2442 798
59 377 2392 436
60 222 2372 388
61 1068 2346 864
62 399 2251 497
63 547 2230 449
64 668 2225 919
65 451 2220 536
66 724 2205 673
67 853 2193 837
68 434 2116 534
69 730 2102 845
70 612 2099 626
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) weekdag zaterdag
200.72174 -0.01949 0.93302
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-1684.3 -227.5 -136.0 135.3 2442.1
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 200.72174 95.96994 2.092 0.0403 *
weekdag -0.01949 0.02081 -0.936 0.3524
zaterdag 0.93302 0.05035 18.532 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 644.3 on 67 degrees of freedom
Multiple R-squared: 0.9807, Adjusted R-squared: 0.9802
F-statistic: 1705 on 2 and 67 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.9999994 1.201596e-06 6.007978e-07
[2,] 1.0000000 6.266405e-12 3.133203e-12
[3,] 1.0000000 3.305516e-13 1.652758e-13
[4,] 1.0000000 1.682398e-12 8.411989e-13
[5,] 1.0000000 3.324329e-14 1.662164e-14
[6,] 1.0000000 1.078429e-13 5.392143e-14
[7,] 1.0000000 2.697733e-13 1.348867e-13
[8,] 1.0000000 9.162529e-13 4.581265e-13
[9,] 1.0000000 3.733569e-12 1.866785e-12
[10,] 1.0000000 3.239384e-12 1.619692e-12
[11,] 1.0000000 4.701308e-12 2.350654e-12
[12,] 1.0000000 8.551020e-12 4.275510e-12
[13,] 1.0000000 7.459107e-12 3.729553e-12
[14,] 1.0000000 8.321189e-12 4.160595e-12
[15,] 1.0000000 2.761699e-11 1.380849e-11
[16,] 1.0000000 1.213752e-12 6.068760e-13
[17,] 1.0000000 2.617730e-12 1.308865e-12
[18,] 1.0000000 4.452667e-13 2.226334e-13
[19,] 1.0000000 1.406360e-12 7.031802e-13
[20,] 1.0000000 4.664994e-12 2.332497e-12
[21,] 1.0000000 1.170146e-11 5.850728e-12
[22,] 1.0000000 3.478378e-11 1.739189e-11
[23,] 1.0000000 1.831525e-11 9.157623e-12
[24,] 1.0000000 3.224161e-11 1.612080e-11
[25,] 1.0000000 1.080571e-10 5.402856e-11
[26,] 1.0000000 3.478421e-10 1.739211e-10
[27,] 1.0000000 1.054839e-09 5.274197e-10
[28,] 1.0000000 2.988786e-09 1.494393e-09
[29,] 1.0000000 8.355484e-09 4.177742e-09
[30,] 1.0000000 2.490620e-08 1.245310e-08
[31,] 1.0000000 7.321986e-08 3.660993e-08
[32,] 0.9999999 1.957287e-07 9.786434e-08
[33,] 0.9999998 4.816294e-07 2.408147e-07
[34,] 0.9999997 6.251020e-07 3.125510e-07
[35,] 0.9999992 1.506265e-06 7.531323e-07
[36,] 0.9999982 3.642598e-06 1.821299e-06
[37,] 0.9999982 3.695042e-06 1.847521e-06
[38,] 0.9999952 9.698830e-06 4.849415e-06
[39,] 0.9999946 1.072480e-05 5.362400e-06
[40,] 0.9999863 2.744848e-05 1.372424e-05
[41,] 0.9999956 8.741993e-06 4.370997e-06
[42,] 0.9999907 1.869847e-05 9.349235e-06
[43,] 0.9999880 2.394235e-05 1.197118e-05
[44,] 0.9999740 5.192225e-05 2.596112e-05
[45,] 0.9999717 5.661702e-05 2.830851e-05
[46,] 0.9999944 1.113708e-05 5.568542e-06
[47,] 0.9999879 2.414294e-05 1.207147e-05
[48,] 0.9999630 7.409105e-05 3.704553e-05
[49,] 0.9999563 8.740510e-05 4.370255e-05
[50,] 0.9999385 1.229367e-04 6.146833e-05
[51,] 0.9998456 3.088112e-04 1.544056e-04
[52,] 0.9995972 8.055821e-04 4.027911e-04
[53,] 0.9999940 1.191932e-05 5.959662e-06
[54,] 0.9999688 6.236670e-05 3.118335e-05
[55,] 0.9999476 1.048761e-04 5.243805e-05
[56,] 0.9999297 1.406087e-04 7.030434e-05
[57,] 0.9996836 6.327275e-04 3.163637e-04
[58,] 0.9983395 3.321056e-03 1.660528e-03
[59,] 0.9980368 3.926317e-03 1.963158e-03
> postscript(file="/var/wessaorg/rcomp/tmp/1cw0n1322138072.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> plot(x[,1], type='l', main='Actuals and Interpolation', ylab='value of Actuals and Interpolation (dots)', xlab='time or index')
> points(x[,1]-mysum$resid)
> grid()
> dev.off()
null device
1
> postscript(file="/var/wessaorg/rcomp/tmp/2gkez1322138072.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> plot(mysum$resid, type='b', pch=19, main='Residuals', ylab='value of Residuals', xlab='time or index')
> grid()
> dev.off()
null device
1
> postscript(file="/var/wessaorg/rcomp/tmp/3tmgu1322138072.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> hist(mysum$resid, main='Residual Histogram', xlab='values of Residuals')
> grid()
> dev.off()
null device
1
> postscript(file="/var/wessaorg/rcomp/tmp/48csd1322138072.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> densityplot(~mysum$resid,col='black',main='Residual Density Plot', xlab='values of Residuals')
> dev.off()
null device
1
> postscript(file="/var/wessaorg/rcomp/tmp/579ca1322138072.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> qqnorm(mysum$resid, main='Residual Normal Q-Q Plot')
> qqline(mysum$resid)
> grid()
> dev.off()
null device
1
> (myerror <- as.ts(mysum$resid))
Time Series:
Start = 1
End = 70
Frequency = 1
1 2 3 4 5 6
-835.25213 -1393.25436 1714.33375 554.68539 -724.46838 -1684.25890
7 8 9 10 11 12
2442.10459 -447.78164 604.89135 2019.48693 520.93208 287.42618
13 14 15 16 17 18
609.51425 626.84029 -197.28558 -171.82897 -129.25138 -502.10586
19 20 21 22 23 24
333.73632 144.77637 -866.32277 -352.89400 1102.37620 -209.28187
25 26 27 28 29 30
-133.08249 -345.81776 -94.73003 837.84517 -336.29975 185.10995
31 32 33 34 35 36
-88.48818 -227.89827 -234.73094 -139.70879 9.38116 -72.14687
37 38 39 40 41 42
-21.50168 -274.92316 -435.68588 -133.25491 -150.80750 460.71847
43 44 45 46 47 48
-11.55320 338.64368 -192.80275 -721.54319 -178.75749 -268.26414
49 50 51 52 53 54
72.65382 -175.39989 765.32942 -226.37933 -25.48180 -225.93098
55 56 57 58 59 60
-138.72378 -154.90554 35.74579 753.31562 -183.90401 -294.50862
61 62 63 64 65 66
106.86513 -221.56636 -29.19045 -346.80931 -206.55845 -61.67510
67 68 69 70
-85.92495 -223.71920 -218.16259 -131.88874
> postscript(file="/var/wessaorg/rcomp/tmp/6i1ya1322138072.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> dum <- cbind(lag(myerror,k=1),myerror)
> dum
Time Series:
Start = 0
End = 70
Frequency = 1
lag(myerror, k = 1) myerror
0 -835.25213 NA
1 -1393.25436 -835.25213
2 1714.33375 -1393.25436
3 554.68539 1714.33375
4 -724.46838 554.68539
5 -1684.25890 -724.46838
6 2442.10459 -1684.25890
7 -447.78164 2442.10459
8 604.89135 -447.78164
9 2019.48693 604.89135
10 520.93208 2019.48693
11 287.42618 520.93208
12 609.51425 287.42618
13 626.84029 609.51425
14 -197.28558 626.84029
15 -171.82897 -197.28558
16 -129.25138 -171.82897
17 -502.10586 -129.25138
18 333.73632 -502.10586
19 144.77637 333.73632
20 -866.32277 144.77637
21 -352.89400 -866.32277
22 1102.37620 -352.89400
23 -209.28187 1102.37620
24 -133.08249 -209.28187
25 -345.81776 -133.08249
26 -94.73003 -345.81776
27 837.84517 -94.73003
28 -336.29975 837.84517
29 185.10995 -336.29975
30 -88.48818 185.10995
31 -227.89827 -88.48818
32 -234.73094 -227.89827
33 -139.70879 -234.73094
34 9.38116 -139.70879
35 -72.14687 9.38116
36 -21.50168 -72.14687
37 -274.92316 -21.50168
38 -435.68588 -274.92316
39 -133.25491 -435.68588
40 -150.80750 -133.25491
41 460.71847 -150.80750
42 -11.55320 460.71847
43 338.64368 -11.55320
44 -192.80275 338.64368
45 -721.54319 -192.80275
46 -178.75749 -721.54319
47 -268.26414 -178.75749
48 72.65382 -268.26414
49 -175.39989 72.65382
50 765.32942 -175.39989
51 -226.37933 765.32942
52 -25.48180 -226.37933
53 -225.93098 -25.48180
54 -138.72378 -225.93098
55 -154.90554 -138.72378
56 35.74579 -154.90554
57 753.31562 35.74579
58 -183.90401 753.31562
59 -294.50862 -183.90401
60 106.86513 -294.50862
61 -221.56636 106.86513
62 -29.19045 -221.56636
63 -346.80931 -29.19045
64 -206.55845 -346.80931
65 -61.67510 -206.55845
66 -85.92495 -61.67510
67 -223.71920 -85.92495
68 -218.16259 -223.71920
69 -131.88874 -218.16259
70 NA -131.88874
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -1393.25436 -835.25213
[2,] 1714.33375 -1393.25436
[3,] 554.68539 1714.33375
[4,] -724.46838 554.68539
[5,] -1684.25890 -724.46838
[6,] 2442.10459 -1684.25890
[7,] -447.78164 2442.10459
[8,] 604.89135 -447.78164
[9,] 2019.48693 604.89135
[10,] 520.93208 2019.48693
[11,] 287.42618 520.93208
[12,] 609.51425 287.42618
[13,] 626.84029 609.51425
[14,] -197.28558 626.84029
[15,] -171.82897 -197.28558
[16,] -129.25138 -171.82897
[17,] -502.10586 -129.25138
[18,] 333.73632 -502.10586
[19,] 144.77637 333.73632
[20,] -866.32277 144.77637
[21,] -352.89400 -866.32277
[22,] 1102.37620 -352.89400
[23,] -209.28187 1102.37620
[24,] -133.08249 -209.28187
[25,] -345.81776 -133.08249
[26,] -94.73003 -345.81776
[27,] 837.84517 -94.73003
[28,] -336.29975 837.84517
[29,] 185.10995 -336.29975
[30,] -88.48818 185.10995
[31,] -227.89827 -88.48818
[32,] -234.73094 -227.89827
[33,] -139.70879 -234.73094
[34,] 9.38116 -139.70879
[35,] -72.14687 9.38116
[36,] -21.50168 -72.14687
[37,] -274.92316 -21.50168
[38,] -435.68588 -274.92316
[39,] -133.25491 -435.68588
[40,] -150.80750 -133.25491
[41,] 460.71847 -150.80750
[42,] -11.55320 460.71847
[43,] 338.64368 -11.55320
[44,] -192.80275 338.64368
[45,] -721.54319 -192.80275
[46,] -178.75749 -721.54319
[47,] -268.26414 -178.75749
[48,] 72.65382 -268.26414
[49,] -175.39989 72.65382
[50,] 765.32942 -175.39989
[51,] -226.37933 765.32942
[52,] -25.48180 -226.37933
[53,] -225.93098 -25.48180
[54,] -138.72378 -225.93098
[55,] -154.90554 -138.72378
[56,] 35.74579 -154.90554
[57,] 753.31562 35.74579
[58,] -183.90401 753.31562
[59,] -294.50862 -183.90401
[60,] 106.86513 -294.50862
[61,] -221.56636 106.86513
[62,] -29.19045 -221.56636
[63,] -346.80931 -29.19045
[64,] -206.55845 -346.80931
[65,] -61.67510 -206.55845
[66,] -85.92495 -61.67510
[67,] -223.71920 -85.92495
[68,] -218.16259 -223.71920
[69,] -131.88874 -218.16259
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -1393.25436 -835.25213
2 1714.33375 -1393.25436
3 554.68539 1714.33375
4 -724.46838 554.68539
5 -1684.25890 -724.46838
6 2442.10459 -1684.25890
7 -447.78164 2442.10459
8 604.89135 -447.78164
9 2019.48693 604.89135
10 520.93208 2019.48693
11 287.42618 520.93208
12 609.51425 287.42618
13 626.84029 609.51425
14 -197.28558 626.84029
15 -171.82897 -197.28558
16 -129.25138 -171.82897
17 -502.10586 -129.25138
18 333.73632 -502.10586
19 144.77637 333.73632
20 -866.32277 144.77637
21 -352.89400 -866.32277
22 1102.37620 -352.89400
23 -209.28187 1102.37620
24 -133.08249 -209.28187
25 -345.81776 -133.08249
26 -94.73003 -345.81776
27 837.84517 -94.73003
28 -336.29975 837.84517
29 185.10995 -336.29975
30 -88.48818 185.10995
31 -227.89827 -88.48818
32 -234.73094 -227.89827
33 -139.70879 -234.73094
34 9.38116 -139.70879
35 -72.14687 9.38116
36 -21.50168 -72.14687
37 -274.92316 -21.50168
38 -435.68588 -274.92316
39 -133.25491 -435.68588
40 -150.80750 -133.25491
41 460.71847 -150.80750
42 -11.55320 460.71847
43 338.64368 -11.55320
44 -192.80275 338.64368
45 -721.54319 -192.80275
46 -178.75749 -721.54319
47 -268.26414 -178.75749
48 72.65382 -268.26414
49 -175.39989 72.65382
50 765.32942 -175.39989
51 -226.37933 765.32942
52 -25.48180 -226.37933
53 -225.93098 -25.48180
54 -138.72378 -225.93098
55 -154.90554 -138.72378
56 35.74579 -154.90554
57 753.31562 35.74579
58 -183.90401 753.31562
59 -294.50862 -183.90401
60 106.86513 -294.50862
61 -221.56636 106.86513
62 -29.19045 -221.56636
63 -346.80931 -29.19045
64 -206.55845 -346.80931
65 -61.67510 -206.55845
66 -85.92495 -61.67510
67 -223.71920 -85.92495
68 -218.16259 -223.71920
69 -131.88874 -218.16259
> 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/wessaorg/rcomp/tmp/7tviv1322138072.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> acf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Autocorrelation Function')
> grid()
> dev.off()
null device
1
> postscript(file="/var/wessaorg/rcomp/tmp/8b96n1322138072.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> pacf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Partial Autocorrelation Function')
> grid()
> dev.off()
null device
1
> postscript(file="/var/wessaorg/rcomp/tmp/957bm1322138072.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> opar <- par(mfrow = c(2,2), oma = c(0, 0, 1.1, 0))
> plot(mylm, las = 1, sub='Residual Diagnostics')
> par(opar)
> dev.off()
null device
1
> if (n > n25) {
+ postscript(file="/var/wessaorg/rcomp/tmp/10rplk1322138072.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
+ plot(kp3:nmkm3,gqarr[,2], main='Goldfeld-Quandt test',ylab='2-sided p-value',xlab='breakpoint')
+ grid()
+ dev.off()
+ }
null device
1
>
> #Note: the /var/wessaorg/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/wessaorg/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/wessaorg/rcomp/tmp/11c2fl1322138072.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/wessaorg/rcomp/tmp/12e8ym1322138072.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/wessaorg/rcomp/tmp/13hbx01322138072.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/wessaorg/rcomp/tmp/14cva51322138072.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/wessaorg/rcomp/tmp/15rs121322138072.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/wessaorg/rcomp/tmp/1627p71322138072.tab")
+ }
>
> try(system("convert tmp/1cw0n1322138072.ps tmp/1cw0n1322138072.png",intern=TRUE))
character(0)
> try(system("convert tmp/2gkez1322138072.ps tmp/2gkez1322138072.png",intern=TRUE))
character(0)
> try(system("convert tmp/3tmgu1322138072.ps tmp/3tmgu1322138072.png",intern=TRUE))
character(0)
> try(system("convert tmp/48csd1322138072.ps tmp/48csd1322138072.png",intern=TRUE))
character(0)
> try(system("convert tmp/579ca1322138072.ps tmp/579ca1322138072.png",intern=TRUE))
character(0)
> try(system("convert tmp/6i1ya1322138072.ps tmp/6i1ya1322138072.png",intern=TRUE))
character(0)
> try(system("convert tmp/7tviv1322138072.ps tmp/7tviv1322138072.png",intern=TRUE))
character(0)
> try(system("convert tmp/8b96n1322138072.ps tmp/8b96n1322138072.png",intern=TRUE))
character(0)
> try(system("convert tmp/957bm1322138072.ps tmp/957bm1322138072.png",intern=TRUE))
character(0)
> try(system("convert tmp/10rplk1322138072.ps tmp/10rplk1322138072.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.247 0.542 3.899