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(210907
+ ,56
+ ,79
+ ,94
+ ,0
+ ,2
+ ,179321
+ ,89
+ ,108
+ ,103
+ ,0
+ ,4
+ ,149061
+ ,44
+ ,43
+ ,93
+ ,0
+ ,0
+ ,237213
+ ,84
+ ,78
+ ,123
+ ,0
+ ,0
+ ,173326
+ ,88
+ ,86
+ ,148
+ ,0
+ ,-4
+ ,133131
+ ,55
+ ,44
+ ,90
+ ,0
+ ,4
+ ,258873
+ ,60
+ ,104
+ ,124
+ ,0
+ ,4
+ ,324799
+ ,154
+ ,158
+ ,168
+ ,0
+ ,0
+ ,230964
+ ,53
+ ,102
+ ,115
+ ,0
+ ,-1
+ ,236785
+ ,119
+ ,77
+ ,71
+ ,0
+ ,0
+ ,344297
+ ,75
+ ,80
+ ,108
+ ,0
+ ,1
+ ,174724
+ ,92
+ ,123
+ ,120
+ ,0
+ ,0
+ ,174415
+ ,100
+ ,73
+ ,114
+ ,0
+ ,3
+ ,223632
+ ,73
+ ,105
+ ,120
+ ,0
+ ,-1
+ ,294424
+ ,77
+ ,107
+ ,124
+ ,0
+ ,4
+ ,325107
+ ,99
+ ,84
+ ,126
+ ,0
+ ,3
+ ,106408
+ ,30
+ ,33
+ ,37
+ ,0
+ ,1
+ ,96560
+ ,76
+ ,42
+ ,38
+ ,1
+ ,0
+ ,265769
+ ,146
+ ,96
+ ,120
+ ,0
+ ,-2
+ ,269651
+ ,67
+ ,106
+ ,93
+ ,0
+ ,-3
+ ,149112
+ ,56
+ ,56
+ ,95
+ ,0
+ ,-4
+ ,152871
+ ,58
+ ,59
+ ,90
+ ,0
+ ,2
+ ,362301
+ ,119
+ ,76
+ ,110
+ ,0
+ ,2
+ ,183167
+ ,66
+ ,91
+ ,138
+ ,0
+ ,-4
+ ,277965
+ ,89
+ ,115
+ ,133
+ ,0
+ ,3
+ ,218946
+ ,41
+ ,76
+ ,96
+ ,0
+ ,2
+ ,244052
+ ,68
+ ,101
+ ,164
+ ,0
+ ,2
+ ,341570
+ ,168
+ ,94
+ ,78
+ ,1
+ ,0
+ ,233328
+ ,132
+ ,92
+ ,102
+ ,0
+ ,5
+ ,206161
+ ,71
+ ,75
+ ,99
+ ,0
+ ,-2
+ ,311473
+ ,112
+ ,128
+ ,129
+ ,0
+ ,0
+ ,207176
+ ,70
+ ,56
+ ,114
+ ,0
+ ,-2
+ ,196553
+ ,57
+ ,41
+ ,99
+ ,0
+ ,-3
+ ,143246
+ ,103
+ ,67
+ ,104
+ ,0
+ ,2
+ ,182192
+ ,52
+ ,77
+ ,138
+ ,0
+ ,2
+ ,194979
+ ,62
+ ,66
+ ,151
+ ,0
+ ,2
+ ,167488
+ ,45
+ ,69
+ ,72
+ ,0
+ ,0
+ ,143756
+ ,46
+ ,105
+ ,120
+ ,0
+ ,4
+ ,275541
+ ,63
+ ,116
+ ,115
+ ,0
+ ,4
+ ,152299
+ ,53
+ ,62
+ ,98
+ ,0
+ ,2
+ ,193339
+ ,78
+ ,100
+ ,71
+ ,0
+ ,2
+ ,130585
+ ,46
+ ,67
+ ,107
+ ,0
+ ,-4
+ ,112611
+ ,41
+ ,46
+ ,73
+ ,1
+ ,3
+ ,148446
+ ,91
+ ,135
+ ,129
+ ,0
+ ,3
+ ,182079
+ ,63
+ ,124
+ ,118
+ ,0
+ ,2
+ ,243060
+ ,63
+ ,58
+ ,104
+ ,0
+ ,-1
+ ,162765
+ ,32
+ ,68
+ ,107
+ ,0
+ ,-3
+ ,85574
+ ,34
+ ,37
+ ,36
+ ,1
+ ,0
+ ,225060
+ ,93
+ ,93
+ ,139
+ ,0
+ ,1
+ ,133328
+ ,55
+ ,56
+ ,56
+ ,1
+ ,-3
+ ,100750
+ ,72
+ ,83
+ ,93
+ ,0
+ ,3
+ ,101523
+ ,42
+ ,59
+ ,87
+ ,1
+ ,0
+ ,243511
+ ,71
+ ,133
+ ,110
+ ,0
+ ,0
+ ,152474
+ ,65
+ ,106
+ ,83
+ ,0
+ ,0
+ ,132487
+ ,41
+ ,71
+ ,98
+ ,0
+ ,3
+ ,317394
+ ,86
+ ,116
+ ,82
+ ,0
+ ,-3
+ ,244749
+ ,95
+ ,98
+ ,115
+ ,0
+ ,0
+ ,184510
+ ,49
+ ,64
+ ,140
+ ,0
+ ,-4
+ ,128423
+ ,64
+ ,32
+ ,120
+ ,0
+ ,2
+ ,97839
+ ,38
+ ,25
+ ,66
+ ,0
+ ,-1
+ ,172494
+ ,52
+ ,46
+ ,139
+ ,0
+ ,3
+ ,229242
+ ,247
+ ,63
+ ,119
+ ,0
+ ,2
+ ,351619
+ ,139
+ ,95
+ ,141
+ ,0
+ ,5
+ ,324598
+ ,110
+ ,113
+ ,133
+ ,0
+ ,2
+ ,195838
+ ,67
+ ,111
+ ,98
+ ,0
+ ,-2
+ ,254488
+ ,83
+ ,120
+ ,117
+ ,0
+ ,0
+ ,199476
+ ,70
+ ,87
+ ,105
+ ,0
+ ,3
+ ,92499
+ ,32
+ ,25
+ ,55
+ ,1
+ ,-2
+ ,224330
+ ,83
+ ,131
+ ,132
+ ,0
+ ,0
+ ,181633
+ ,70
+ ,47
+ ,73
+ ,0
+ ,6
+ ,271856
+ ,103
+ ,109
+ ,86
+ ,0
+ ,-3
+ ,95227
+ ,34
+ ,37
+ ,48
+ ,0
+ ,3
+ ,98146
+ ,40
+ ,15
+ ,48
+ ,1
+ ,0
+ ,118612
+ ,46
+ ,54
+ ,43
+ ,1
+ ,-2
+ ,65475
+ ,18
+ ,16
+ ,46
+ ,1
+ ,1
+ ,108446
+ ,60
+ ,22
+ ,65
+ ,1
+ ,0
+ ,121848
+ ,39
+ ,37
+ ,52
+ ,1
+ ,2
+ ,76302
+ ,31
+ ,29
+ ,68
+ ,1
+ ,2
+ ,98104
+ ,54
+ ,55
+ ,47
+ ,1
+ ,-3
+ ,30989
+ ,14
+ ,5
+ ,41
+ ,1
+ ,-2
+ ,31774
+ ,23
+ ,0
+ ,47
+ ,1
+ ,1
+ ,150580
+ ,77
+ ,27
+ ,71
+ ,1
+ ,-4
+ ,54157
+ ,19
+ ,37
+ ,30
+ ,1
+ ,0
+ ,59382
+ ,49
+ ,29
+ ,24
+ ,1
+ ,1
+ ,84105
+ ,20
+ ,17
+ ,63
+ ,1
+ ,0)
+ ,dim=c(6
+ ,85)
+ ,dimnames=list(c('time'
+ ,'logins'
+ ,'BC'
+ ,'LFM'
+ ,'Course'
+ ,'Totaal')
+ ,1:85))
> y <- array(NA,dim=c(6,85),dimnames=list(c('time','logins','BC','LFM','Course','Totaal'),1:85))
> 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
> 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
time logins BC LFM Course Totaal
1 210907 56 79 94 0 2
2 179321 89 108 103 0 4
3 149061 44 43 93 0 0
4 237213 84 78 123 0 0
5 173326 88 86 148 0 -4
6 133131 55 44 90 0 4
7 258873 60 104 124 0 4
8 324799 154 158 168 0 0
9 230964 53 102 115 0 -1
10 236785 119 77 71 0 0
11 344297 75 80 108 0 1
12 174724 92 123 120 0 0
13 174415 100 73 114 0 3
14 223632 73 105 120 0 -1
15 294424 77 107 124 0 4
16 325107 99 84 126 0 3
17 106408 30 33 37 0 1
18 96560 76 42 38 1 0
19 265769 146 96 120 0 -2
20 269651 67 106 93 0 -3
21 149112 56 56 95 0 -4
22 152871 58 59 90 0 2
23 362301 119 76 110 0 2
24 183167 66 91 138 0 -4
25 277965 89 115 133 0 3
26 218946 41 76 96 0 2
27 244052 68 101 164 0 2
28 341570 168 94 78 1 0
29 233328 132 92 102 0 5
30 206161 71 75 99 0 -2
31 311473 112 128 129 0 0
32 207176 70 56 114 0 -2
33 196553 57 41 99 0 -3
34 143246 103 67 104 0 2
35 182192 52 77 138 0 2
36 194979 62 66 151 0 2
37 167488 45 69 72 0 0
38 143756 46 105 120 0 4
39 275541 63 116 115 0 4
40 152299 53 62 98 0 2
41 193339 78 100 71 0 2
42 130585 46 67 107 0 -4
43 112611 41 46 73 1 3
44 148446 91 135 129 0 3
45 182079 63 124 118 0 2
46 243060 63 58 104 0 -1
47 162765 32 68 107 0 -3
48 85574 34 37 36 1 0
49 225060 93 93 139 0 1
50 133328 55 56 56 1 -3
51 100750 72 83 93 0 3
52 101523 42 59 87 1 0
53 243511 71 133 110 0 0
54 152474 65 106 83 0 0
55 132487 41 71 98 0 3
56 317394 86 116 82 0 -3
57 244749 95 98 115 0 0
58 184510 49 64 140 0 -4
59 128423 64 32 120 0 2
60 97839 38 25 66 0 -1
61 172494 52 46 139 0 3
62 229242 247 63 119 0 2
63 351619 139 95 141 0 5
64 324598 110 113 133 0 2
65 195838 67 111 98 0 -2
66 254488 83 120 117 0 0
67 199476 70 87 105 0 3
68 92499 32 25 55 1 -2
69 224330 83 131 132 0 0
70 181633 70 47 73 0 6
71 271856 103 109 86 0 -3
72 95227 34 37 48 0 3
73 98146 40 15 48 1 0
74 118612 46 54 43 1 -2
75 65475 18 16 46 1 1
76 108446 60 22 65 1 0
77 121848 39 37 52 1 2
78 76302 31 29 68 1 2
79 98104 54 55 47 1 -3
80 30989 14 5 41 1 -2
81 31774 23 0 47 1 1
82 150580 77 27 71 1 -4
83 54157 19 37 30 1 0
84 59382 49 29 24 1 1
85 84105 20 17 63 1 0
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) logins BC LFM Course Totaal
38756.5 818.9 796.8 376.5 -22382.5 -743.9
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-118740 -24618 -781 22061 140456
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 38756.5 25594.4 1.514 0.133952
logins 818.9 167.2 4.898 5.04e-06 ***
BC 796.8 218.9 3.639 0.000486 ***
LFM 376.5 247.6 1.521 0.132282
Course -22382.5 17971.1 -1.245 0.216641
Totaal -743.9 2127.6 -0.350 0.727544
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 46990 on 79 degrees of freedom
Multiple R-squared: 0.6744, Adjusted R-squared: 0.6538
F-statistic: 32.72 on 5 and 79 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.5809897 8.380206e-01 4.190103e-01
[2,] 0.5474391 9.051219e-01 4.525609e-01
[3,] 0.9664984 6.700318e-02 3.350159e-02
[4,] 0.9832079 3.358413e-02 1.679207e-02
[5,] 0.9774350 4.512994e-02 2.256497e-02
[6,] 0.9602558 7.948837e-02 3.974419e-02
[7,] 0.9621088 7.578232e-02 3.789116e-02
[8,] 0.9839842 3.203150e-02 1.601575e-02
[9,] 0.9737979 5.240427e-02 2.620213e-02
[10,] 0.9596246 8.075078e-02 4.037539e-02
[11,] 0.9407992 1.184015e-01 5.920077e-02
[12,] 0.9501539 9.969218e-02 4.984609e-02
[13,] 0.9294479 1.411043e-01 7.055214e-02
[14,] 0.9079565 1.840870e-01 9.204348e-02
[15,] 0.9866349 2.673012e-02 1.336506e-02
[16,] 0.9824675 3.506508e-02 1.753254e-02
[17,] 0.9750933 4.981349e-02 2.490675e-02
[18,] 0.9743029 5.139417e-02 2.569708e-02
[19,] 0.9621302 7.573954e-02 3.786977e-02
[20,] 0.9830795 3.384093e-02 1.692047e-02
[21,] 0.9823293 3.534146e-02 1.767073e-02
[22,] 0.9738489 5.230220e-02 2.615110e-02
[23,] 0.9680197 6.396056e-02 3.198028e-02
[24,] 0.9568871 8.622582e-02 4.311291e-02
[25,] 0.9506484 9.870327e-02 4.935164e-02
[26,] 0.9699974 6.000512e-02 3.000256e-02
[27,] 0.9585112 8.297757e-02 4.148879e-02
[28,] 0.9417870 1.164260e-01 5.821300e-02
[29,] 0.9213727 1.572545e-01 7.862726e-02
[30,] 0.9348158 1.303684e-01 6.518422e-02
[31,] 0.9433607 1.132787e-01 5.663933e-02
[32,] 0.9245993 1.508014e-01 7.540071e-02
[33,] 0.9025384 1.949233e-01 9.746163e-02
[34,] 0.8997093 2.005814e-01 1.002907e-01
[35,] 0.8711727 2.576545e-01 1.288273e-01
[36,] 0.9727927 5.441466e-02 2.720733e-02
[37,] 0.9743598 5.128036e-02 2.564018e-02
[38,] 0.9853059 2.938825e-02 1.469412e-02
[39,] 0.9776306 4.473883e-02 2.236942e-02
[40,] 0.9669322 6.613554e-02 3.306777e-02
[41,] 0.9543744 9.125119e-02 4.562560e-02
[42,] 0.9350605 1.298791e-01 6.493954e-02
[43,] 0.9828755 3.424892e-02 1.712446e-02
[44,] 0.9818220 3.635608e-02 1.817804e-02
[45,] 0.9744297 5.114061e-02 2.557030e-02
[46,] 0.9836760 3.264801e-02 1.632400e-02
[47,] 0.9834170 3.316604e-02 1.658302e-02
[48,] 0.9969486 6.102783e-03 3.051391e-03
[49,] 0.9946058 1.078841e-02 5.394203e-03
[50,] 0.9906646 1.867085e-02 9.335425e-03
[51,] 0.9884110 2.317801e-02 1.158901e-02
[52,] 0.9814670 3.706596e-02 1.853298e-02
[53,] 0.9740997 5.180057e-02 2.590028e-02
[54,] 0.9999546 9.088765e-05 4.544383e-05
[55,] 0.9999125 1.750453e-04 8.752264e-05
[56,] 0.9999301 1.397846e-04 6.989228e-05
[57,] 0.9998305 3.389376e-04 1.694688e-04
[58,] 0.9996307 7.386935e-04 3.693467e-04
[59,] 0.9990340 1.932029e-03 9.660147e-04
[60,] 0.9979025 4.194921e-03 2.097460e-03
[61,] 0.9995610 8.779812e-04 4.389906e-04
[62,] 0.9987794 2.441150e-03 1.220575e-03
[63,] 0.9967243 6.551347e-03 3.275673e-03
[64,] 0.9909869 1.802626e-02 9.013131e-03
[65,] 0.9869072 2.618567e-02 1.309283e-02
[66,] 0.9703748 5.925033e-02 2.962517e-02
[67,] 0.9355264 1.289471e-01 6.447357e-02
[68,] 0.8484896 3.030207e-01 1.515104e-01
> postscript(file="/var/wessaorg/rcomp/tmp/1n2vg1324636345.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/21aqc1324636345.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/3issx1324636345.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/42qli1324636345.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/5hwp61324636345.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 = 85
Frequency = 1
1 2 3 4 5 6
29437.8096 -54180.7221 4991.8789 21203.0836 -64722.8459 -16638.0479
7 8 9 10 11 12
44400.2157 -29222.9032 23486.7512 12489.0932 140455.8728 -82562.8789
13 14 15 16 17 18
-45093.4926 -4496.9645 63638.9896 93134.4154 3601.7271 -29826.1876
19 20 21 22 23 24
-15715.0424 54317.6076 -18871.0697 -12794.3112 125604.6811 -37083.8146
25 26 27 28 29 30
26845.8207 51398.0922 8868.8296 83347.7921 -21519.0081 10736.7052
31 32 33 34 35 36
30434.6964 22061.4003 38940.4265 -70917.1241 -10975.5312 -2508.2807
37 38 39 40 41 42
9790.9475 -58542.3325 52438.8879 -14674.2766 -14218.7717 -42491.7513
43 44 45 46 47 48
753.2365 -118740.4986 -50014.7239 66593.6344 1100.4317 -1679.9650
49 50 51 52 53 54
-15552.7419 3975.2061 -95888.8446 -29014.9994 -780.6819 -55224.5105
55 56 57 58 59 60
-31086.2022 82674.9624 6807.5029 -1058.9439 -31938.9764 -17551.8495
61 62 63 64 65 66
4393.9905 -105308.6286 73964.1581 57130.8762 -24618.0983 8091.4831
67 68 69 70 71 72
-3230.5424 7802.2955 -36479.1551 25078.4884 27286.3870 -16696.2800
73 74 75 76 77 78
18989.1403 3862.0408 5034.8737 931.8578 25962.5271 -12682.3532
79 80 81 82 83 84
-26244.2442 -17759.7631 -20388.8566 19925.3276 -18553.7243 -28519.4155
85
14085.2115
> postscript(file="/var/wessaorg/rcomp/tmp/61vgv1324636345.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 = 85
Frequency = 1
lag(myerror, k = 1) myerror
0 29437.8096 NA
1 -54180.7221 29437.8096
2 4991.8789 -54180.7221
3 21203.0836 4991.8789
4 -64722.8459 21203.0836
5 -16638.0479 -64722.8459
6 44400.2157 -16638.0479
7 -29222.9032 44400.2157
8 23486.7512 -29222.9032
9 12489.0932 23486.7512
10 140455.8728 12489.0932
11 -82562.8789 140455.8728
12 -45093.4926 -82562.8789
13 -4496.9645 -45093.4926
14 63638.9896 -4496.9645
15 93134.4154 63638.9896
16 3601.7271 93134.4154
17 -29826.1876 3601.7271
18 -15715.0424 -29826.1876
19 54317.6076 -15715.0424
20 -18871.0697 54317.6076
21 -12794.3112 -18871.0697
22 125604.6811 -12794.3112
23 -37083.8146 125604.6811
24 26845.8207 -37083.8146
25 51398.0922 26845.8207
26 8868.8296 51398.0922
27 83347.7921 8868.8296
28 -21519.0081 83347.7921
29 10736.7052 -21519.0081
30 30434.6964 10736.7052
31 22061.4003 30434.6964
32 38940.4265 22061.4003
33 -70917.1241 38940.4265
34 -10975.5312 -70917.1241
35 -2508.2807 -10975.5312
36 9790.9475 -2508.2807
37 -58542.3325 9790.9475
38 52438.8879 -58542.3325
39 -14674.2766 52438.8879
40 -14218.7717 -14674.2766
41 -42491.7513 -14218.7717
42 753.2365 -42491.7513
43 -118740.4986 753.2365
44 -50014.7239 -118740.4986
45 66593.6344 -50014.7239
46 1100.4317 66593.6344
47 -1679.9650 1100.4317
48 -15552.7419 -1679.9650
49 3975.2061 -15552.7419
50 -95888.8446 3975.2061
51 -29014.9994 -95888.8446
52 -780.6819 -29014.9994
53 -55224.5105 -780.6819
54 -31086.2022 -55224.5105
55 82674.9624 -31086.2022
56 6807.5029 82674.9624
57 -1058.9439 6807.5029
58 -31938.9764 -1058.9439
59 -17551.8495 -31938.9764
60 4393.9905 -17551.8495
61 -105308.6286 4393.9905
62 73964.1581 -105308.6286
63 57130.8762 73964.1581
64 -24618.0983 57130.8762
65 8091.4831 -24618.0983
66 -3230.5424 8091.4831
67 7802.2955 -3230.5424
68 -36479.1551 7802.2955
69 25078.4884 -36479.1551
70 27286.3870 25078.4884
71 -16696.2800 27286.3870
72 18989.1403 -16696.2800
73 3862.0408 18989.1403
74 5034.8737 3862.0408
75 931.8578 5034.8737
76 25962.5271 931.8578
77 -12682.3532 25962.5271
78 -26244.2442 -12682.3532
79 -17759.7631 -26244.2442
80 -20388.8566 -17759.7631
81 19925.3276 -20388.8566
82 -18553.7243 19925.3276
83 -28519.4155 -18553.7243
84 14085.2115 -28519.4155
85 NA 14085.2115
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -54180.7221 29437.8096
[2,] 4991.8789 -54180.7221
[3,] 21203.0836 4991.8789
[4,] -64722.8459 21203.0836
[5,] -16638.0479 -64722.8459
[6,] 44400.2157 -16638.0479
[7,] -29222.9032 44400.2157
[8,] 23486.7512 -29222.9032
[9,] 12489.0932 23486.7512
[10,] 140455.8728 12489.0932
[11,] -82562.8789 140455.8728
[12,] -45093.4926 -82562.8789
[13,] -4496.9645 -45093.4926
[14,] 63638.9896 -4496.9645
[15,] 93134.4154 63638.9896
[16,] 3601.7271 93134.4154
[17,] -29826.1876 3601.7271
[18,] -15715.0424 -29826.1876
[19,] 54317.6076 -15715.0424
[20,] -18871.0697 54317.6076
[21,] -12794.3112 -18871.0697
[22,] 125604.6811 -12794.3112
[23,] -37083.8146 125604.6811
[24,] 26845.8207 -37083.8146
[25,] 51398.0922 26845.8207
[26,] 8868.8296 51398.0922
[27,] 83347.7921 8868.8296
[28,] -21519.0081 83347.7921
[29,] 10736.7052 -21519.0081
[30,] 30434.6964 10736.7052
[31,] 22061.4003 30434.6964
[32,] 38940.4265 22061.4003
[33,] -70917.1241 38940.4265
[34,] -10975.5312 -70917.1241
[35,] -2508.2807 -10975.5312
[36,] 9790.9475 -2508.2807
[37,] -58542.3325 9790.9475
[38,] 52438.8879 -58542.3325
[39,] -14674.2766 52438.8879
[40,] -14218.7717 -14674.2766
[41,] -42491.7513 -14218.7717
[42,] 753.2365 -42491.7513
[43,] -118740.4986 753.2365
[44,] -50014.7239 -118740.4986
[45,] 66593.6344 -50014.7239
[46,] 1100.4317 66593.6344
[47,] -1679.9650 1100.4317
[48,] -15552.7419 -1679.9650
[49,] 3975.2061 -15552.7419
[50,] -95888.8446 3975.2061
[51,] -29014.9994 -95888.8446
[52,] -780.6819 -29014.9994
[53,] -55224.5105 -780.6819
[54,] -31086.2022 -55224.5105
[55,] 82674.9624 -31086.2022
[56,] 6807.5029 82674.9624
[57,] -1058.9439 6807.5029
[58,] -31938.9764 -1058.9439
[59,] -17551.8495 -31938.9764
[60,] 4393.9905 -17551.8495
[61,] -105308.6286 4393.9905
[62,] 73964.1581 -105308.6286
[63,] 57130.8762 73964.1581
[64,] -24618.0983 57130.8762
[65,] 8091.4831 -24618.0983
[66,] -3230.5424 8091.4831
[67,] 7802.2955 -3230.5424
[68,] -36479.1551 7802.2955
[69,] 25078.4884 -36479.1551
[70,] 27286.3870 25078.4884
[71,] -16696.2800 27286.3870
[72,] 18989.1403 -16696.2800
[73,] 3862.0408 18989.1403
[74,] 5034.8737 3862.0408
[75,] 931.8578 5034.8737
[76,] 25962.5271 931.8578
[77,] -12682.3532 25962.5271
[78,] -26244.2442 -12682.3532
[79,] -17759.7631 -26244.2442
[80,] -20388.8566 -17759.7631
[81,] 19925.3276 -20388.8566
[82,] -18553.7243 19925.3276
[83,] -28519.4155 -18553.7243
[84,] 14085.2115 -28519.4155
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -54180.7221 29437.8096
2 4991.8789 -54180.7221
3 21203.0836 4991.8789
4 -64722.8459 21203.0836
5 -16638.0479 -64722.8459
6 44400.2157 -16638.0479
7 -29222.9032 44400.2157
8 23486.7512 -29222.9032
9 12489.0932 23486.7512
10 140455.8728 12489.0932
11 -82562.8789 140455.8728
12 -45093.4926 -82562.8789
13 -4496.9645 -45093.4926
14 63638.9896 -4496.9645
15 93134.4154 63638.9896
16 3601.7271 93134.4154
17 -29826.1876 3601.7271
18 -15715.0424 -29826.1876
19 54317.6076 -15715.0424
20 -18871.0697 54317.6076
21 -12794.3112 -18871.0697
22 125604.6811 -12794.3112
23 -37083.8146 125604.6811
24 26845.8207 -37083.8146
25 51398.0922 26845.8207
26 8868.8296 51398.0922
27 83347.7921 8868.8296
28 -21519.0081 83347.7921
29 10736.7052 -21519.0081
30 30434.6964 10736.7052
31 22061.4003 30434.6964
32 38940.4265 22061.4003
33 -70917.1241 38940.4265
34 -10975.5312 -70917.1241
35 -2508.2807 -10975.5312
36 9790.9475 -2508.2807
37 -58542.3325 9790.9475
38 52438.8879 -58542.3325
39 -14674.2766 52438.8879
40 -14218.7717 -14674.2766
41 -42491.7513 -14218.7717
42 753.2365 -42491.7513
43 -118740.4986 753.2365
44 -50014.7239 -118740.4986
45 66593.6344 -50014.7239
46 1100.4317 66593.6344
47 -1679.9650 1100.4317
48 -15552.7419 -1679.9650
49 3975.2061 -15552.7419
50 -95888.8446 3975.2061
51 -29014.9994 -95888.8446
52 -780.6819 -29014.9994
53 -55224.5105 -780.6819
54 -31086.2022 -55224.5105
55 82674.9624 -31086.2022
56 6807.5029 82674.9624
57 -1058.9439 6807.5029
58 -31938.9764 -1058.9439
59 -17551.8495 -31938.9764
60 4393.9905 -17551.8495
61 -105308.6286 4393.9905
62 73964.1581 -105308.6286
63 57130.8762 73964.1581
64 -24618.0983 57130.8762
65 8091.4831 -24618.0983
66 -3230.5424 8091.4831
67 7802.2955 -3230.5424
68 -36479.1551 7802.2955
69 25078.4884 -36479.1551
70 27286.3870 25078.4884
71 -16696.2800 27286.3870
72 18989.1403 -16696.2800
73 3862.0408 18989.1403
74 5034.8737 3862.0408
75 931.8578 5034.8737
76 25962.5271 931.8578
77 -12682.3532 25962.5271
78 -26244.2442 -12682.3532
79 -17759.7631 -26244.2442
80 -20388.8566 -17759.7631
81 19925.3276 -20388.8566
82 -18553.7243 19925.3276
83 -28519.4155 -18553.7243
84 14085.2115 -28519.4155
> 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/7ry381324636345.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/87epv1324636345.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/9glbi1324636345.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/10v2ve1324636345.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/11vocx1324636345.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/12eyfw1324636345.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/13u52v1324636345.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/140jrs1324636345.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/15iv0q1324636345.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/169qlp1324636345.tab")
+ }
>
> try(system("convert tmp/1n2vg1324636345.ps tmp/1n2vg1324636345.png",intern=TRUE))
character(0)
> try(system("convert tmp/21aqc1324636345.ps tmp/21aqc1324636345.png",intern=TRUE))
character(0)
> try(system("convert tmp/3issx1324636345.ps tmp/3issx1324636345.png",intern=TRUE))
character(0)
> try(system("convert tmp/42qli1324636345.ps tmp/42qli1324636345.png",intern=TRUE))
character(0)
> try(system("convert tmp/5hwp61324636345.ps tmp/5hwp61324636345.png",intern=TRUE))
character(0)
> try(system("convert tmp/61vgv1324636345.ps tmp/61vgv1324636345.png",intern=TRUE))
character(0)
> try(system("convert tmp/7ry381324636345.ps tmp/7ry381324636345.png",intern=TRUE))
character(0)
> try(system("convert tmp/87epv1324636345.ps tmp/87epv1324636345.png",intern=TRUE))
character(0)
> try(system("convert tmp/9glbi1324636345.ps tmp/9glbi1324636345.png",intern=TRUE))
character(0)
> try(system("convert tmp/10v2ve1324636345.ps tmp/10v2ve1324636345.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.538 0.575 4.162