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(180144
+ ,40.6
+ ,173666
+ ,63.6
+ ,165688
+ ,66.8
+ ,161570
+ ,71.5
+ ,156145
+ ,99.4
+ ,153730
+ ,78.2
+ ,182698
+ ,57.2
+ ,200765
+ ,86.5
+ ,176512
+ ,66.1
+ ,166618
+ ,75
+ ,158644
+ ,55
+ ,159585
+ ,66.8
+ ,163095
+ ,41.4
+ ,159044
+ ,53.3
+ ,155511
+ ,71.4
+ ,153745
+ ,68.2
+ ,150569
+ ,84.1
+ ,150605
+ ,94
+ ,179612
+ ,91.4
+ ,194690
+ ,79.9
+ ,189917
+ ,40.7
+ ,184128
+ ,60.3
+ ,175335
+ ,49.1
+ ,179566
+ ,42
+ ,181140
+ ,54.3
+ ,177876
+ ,39.3
+ ,175041
+ ,47.8
+ ,169292
+ ,74.5
+ ,166070
+ ,78.8
+ ,166972
+ ,81.4
+ ,206348
+ ,66
+ ,215706
+ ,88.8
+ ,202108
+ ,54.4
+ ,195411
+ ,75.8
+ ,193111
+ ,51.6
+ ,195198
+ ,53
+ ,198770
+ ,62.7
+ ,194163
+ ,52.3
+ ,190420
+ ,30.5
+ ,189733
+ ,49.9
+ ,186029
+ ,53.8
+ ,191531
+ ,65.3
+ ,232571
+ ,62.7
+ ,243477
+ ,55.4
+ ,227247
+ ,66.2
+ ,217859
+ ,67.2
+ ,208679
+ ,42.4
+ ,213188
+ ,56.3
+ ,216234
+ ,44.9
+ ,213586
+ ,30
+ ,209465
+ ,54.4
+ ,204045
+ ,47.8
+ ,200237
+ ,63.6
+ ,203666
+ ,72.5
+ ,241476
+ ,82.2
+ ,260307
+ ,67.9
+ ,243324
+ ,67.8
+ ,244460
+ ,65.6
+ ,233575
+ ,78.1
+ ,237217
+ ,41.6
+ ,235243
+ ,64.3
+ ,230354
+ ,55.9
+ ,227184
+ ,78.3
+ ,221678
+ ,69.8
+ ,217142
+ ,59.3
+ ,219452
+ ,103.6
+ ,256446
+ ,109.7
+ ,265845
+ ,76.3
+ ,248624
+ ,81.8
+ ,241114
+ ,99.6
+ ,229245
+ ,100.6
+ ,231805
+ ,79.9
+ ,219277
+ ,49.3
+ ,219313
+ ,62.7
+ ,212610
+ ,101.3
+ ,214771
+ ,101.2
+ ,211142
+ ,83.3
+ ,211457
+ ,127.8
+ ,240048
+ ,103.7
+ ,240636
+ ,91.5
+ ,230580
+ ,95.1
+ ,208795
+ ,109
+ ,197922
+ ,132.6
+ ,194596
+ ,79.5)
+ ,dim=c(2
+ ,84)
+ ,dimnames=list(c('werkloosheid'
+ ,'Oceanië')
+ ,1:84))
> y <- array(NA,dim=c(2,84),dimnames=list(c('werkloosheid','Oceanië'),1:84))
> 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
werkloosheid Oceani\353 M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 180144 40.6 1 0 0 0 0 0 0 0 0 0 0 1
2 173666 63.6 0 1 0 0 0 0 0 0 0 0 0 2
3 165688 66.8 0 0 1 0 0 0 0 0 0 0 0 3
4 161570 71.5 0 0 0 1 0 0 0 0 0 0 0 4
5 156145 99.4 0 0 0 0 1 0 0 0 0 0 0 5
6 153730 78.2 0 0 0 0 0 1 0 0 0 0 0 6
7 182698 57.2 0 0 0 0 0 0 1 0 0 0 0 7
8 200765 86.5 0 0 0 0 0 0 0 1 0 0 0 8
9 176512 66.1 0 0 0 0 0 0 0 0 1 0 0 9
10 166618 75.0 0 0 0 0 0 0 0 0 0 1 0 10
11 158644 55.0 0 0 0 0 0 0 0 0 0 0 1 11
12 159585 66.8 0 0 0 0 0 0 0 0 0 0 0 12
13 163095 41.4 1 0 0 0 0 0 0 0 0 0 0 13
14 159044 53.3 0 1 0 0 0 0 0 0 0 0 0 14
15 155511 71.4 0 0 1 0 0 0 0 0 0 0 0 15
16 153745 68.2 0 0 0 1 0 0 0 0 0 0 0 16
17 150569 84.1 0 0 0 0 1 0 0 0 0 0 0 17
18 150605 94.0 0 0 0 0 0 1 0 0 0 0 0 18
19 179612 91.4 0 0 0 0 0 0 1 0 0 0 0 19
20 194690 79.9 0 0 0 0 0 0 0 1 0 0 0 20
21 189917 40.7 0 0 0 0 0 0 0 0 1 0 0 21
22 184128 60.3 0 0 0 0 0 0 0 0 0 1 0 22
23 175335 49.1 0 0 0 0 0 0 0 0 0 0 1 23
24 179566 42.0 0 0 0 0 0 0 0 0 0 0 0 24
25 181140 54.3 1 0 0 0 0 0 0 0 0 0 0 25
26 177876 39.3 0 1 0 0 0 0 0 0 0 0 0 26
27 175041 47.8 0 0 1 0 0 0 0 0 0 0 0 27
28 169292 74.5 0 0 0 1 0 0 0 0 0 0 0 28
29 166070 78.8 0 0 0 0 1 0 0 0 0 0 0 29
30 166972 81.4 0 0 0 0 0 1 0 0 0 0 0 30
31 206348 66.0 0 0 0 0 0 0 1 0 0 0 0 31
32 215706 88.8 0 0 0 0 0 0 0 1 0 0 0 32
33 202108 54.4 0 0 0 0 0 0 0 0 1 0 0 33
34 195411 75.8 0 0 0 0 0 0 0 0 0 1 0 34
35 193111 51.6 0 0 0 0 0 0 0 0 0 0 1 35
36 195198 53.0 0 0 0 0 0 0 0 0 0 0 0 36
37 198770 62.7 1 0 0 0 0 0 0 0 0 0 0 37
38 194163 52.3 0 1 0 0 0 0 0 0 0 0 0 38
39 190420 30.5 0 0 1 0 0 0 0 0 0 0 0 39
40 189733 49.9 0 0 0 1 0 0 0 0 0 0 0 40
41 186029 53.8 0 0 0 0 1 0 0 0 0 0 0 41
42 191531 65.3 0 0 0 0 0 1 0 0 0 0 0 42
43 232571 62.7 0 0 0 0 0 0 1 0 0 0 0 43
44 243477 55.4 0 0 0 0 0 0 0 1 0 0 0 44
45 227247 66.2 0 0 0 0 0 0 0 0 1 0 0 45
46 217859 67.2 0 0 0 0 0 0 0 0 0 1 0 46
47 208679 42.4 0 0 0 0 0 0 0 0 0 0 1 47
48 213188 56.3 0 0 0 0 0 0 0 0 0 0 0 48
49 216234 44.9 1 0 0 0 0 0 0 0 0 0 0 49
50 213586 30.0 0 1 0 0 0 0 0 0 0 0 0 50
51 209465 54.4 0 0 1 0 0 0 0 0 0 0 0 51
52 204045 47.8 0 0 0 1 0 0 0 0 0 0 0 52
53 200237 63.6 0 0 0 0 1 0 0 0 0 0 0 53
54 203666 72.5 0 0 0 0 0 1 0 0 0 0 0 54
55 241476 82.2 0 0 0 0 0 0 1 0 0 0 0 55
56 260307 67.9 0 0 0 0 0 0 0 1 0 0 0 56
57 243324 67.8 0 0 0 0 0 0 0 0 1 0 0 57
58 244460 65.6 0 0 0 0 0 0 0 0 0 1 0 58
59 233575 78.1 0 0 0 0 0 0 0 0 0 0 1 59
60 237217 41.6 0 0 0 0 0 0 0 0 0 0 0 60
61 235243 64.3 1 0 0 0 0 0 0 0 0 0 0 61
62 230354 55.9 0 1 0 0 0 0 0 0 0 0 0 62
63 227184 78.3 0 0 1 0 0 0 0 0 0 0 0 63
64 221678 69.8 0 0 0 1 0 0 0 0 0 0 0 64
65 217142 59.3 0 0 0 0 1 0 0 0 0 0 0 65
66 219452 103.6 0 0 0 0 0 1 0 0 0 0 0 66
67 256446 109.7 0 0 0 0 0 0 1 0 0 0 0 67
68 265845 76.3 0 0 0 0 0 0 0 1 0 0 0 68
69 248624 81.8 0 0 0 0 0 0 0 0 1 0 0 69
70 241114 99.6 0 0 0 0 0 0 0 0 0 1 0 70
71 229245 100.6 0 0 0 0 0 0 0 0 0 0 1 71
72 231805 79.9 0 0 0 0 0 0 0 0 0 0 0 72
73 219277 49.3 1 0 0 0 0 0 0 0 0 0 0 73
74 219313 62.7 0 1 0 0 0 0 0 0 0 0 0 74
75 212610 101.3 0 0 1 0 0 0 0 0 0 0 0 75
76 214771 101.2 0 0 0 1 0 0 0 0 0 0 0 76
77 211142 83.3 0 0 0 0 1 0 0 0 0 0 0 77
78 211457 127.8 0 0 0 0 0 1 0 0 0 0 0 78
79 240048 103.7 0 0 0 0 0 0 1 0 0 0 0 79
80 240636 91.5 0 0 0 0 0 0 0 1 0 0 0 80
81 230580 95.1 0 0 0 0 0 0 0 0 1 0 0 81
82 208795 109.0 0 0 0 0 0 0 0 0 0 1 0 82
83 197922 132.6 0 0 0 0 0 0 0 0 0 0 1 83
84 194596 79.5 0 0 0 0 0 0 0 0 0 0 0 84
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) `Oceani\353` M1 M2 M3
171222.0 -316.7 6054.1 1308.0 -77.0
M4 M5 M6 M7 M8
-2650.9 -5824.6 -865.2 30390.0 39905.4
M9 M10 M11 t
20789.7 14838.7 3021.6 1027.8
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-37781.6 -8103.0 -845.5 9564.3 23427.5
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 171222.03 6714.87 25.499 < 2e-16 ***
`Oceani\353` -316.73 86.11 -3.678 0.000458 ***
M1 6054.05 6795.38 0.891 0.376031
M2 1308.01 6792.54 0.193 0.847857
M3 -77.00 6800.00 -0.011 0.990997
M4 -2650.95 6840.48 -0.388 0.699534
M5 -5824.59 6918.86 -0.842 0.402744
M6 -865.20 7268.03 -0.119 0.905583
M7 30390.02 7055.00 4.308 5.28e-05 ***
M8 39905.36 6959.11 5.734 2.29e-07 ***
M9 20789.65 6792.67 3.061 0.003131 **
M10 14838.71 6961.41 2.132 0.036555 *
M11 3021.62 6846.20 0.441 0.660315
t 1027.80 63.67 16.143 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 12630 on 70 degrees of freedom
Multiple R-squared: 0.8449, Adjusted R-squared: 0.8161
F-statistic: 29.34 on 13 and 70 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.027939397 0.055878795 0.9720606
[2,] 0.017400401 0.034800803 0.9825996
[3,] 0.004811486 0.009622973 0.9951885
[4,] 0.001639075 0.003278149 0.9983609
[5,] 0.027222499 0.054444999 0.9727775
[6,] 0.060223369 0.120446738 0.9397766
[7,] 0.083075531 0.166151062 0.9169245
[8,] 0.076295736 0.152591473 0.9237043
[9,] 0.075854067 0.151708134 0.9241459
[10,] 0.049892046 0.099784092 0.9501080
[11,] 0.033788568 0.067577136 0.9662114
[12,] 0.025595574 0.051191148 0.9744044
[13,] 0.015405500 0.030810999 0.9845945
[14,] 0.012126047 0.024252094 0.9878740
[15,] 0.018350756 0.036701512 0.9816492
[16,] 0.015910126 0.031820252 0.9840899
[17,] 0.015324090 0.030648179 0.9846759
[18,] 0.014821948 0.029643895 0.9851781
[19,] 0.017277738 0.034555476 0.9827223
[20,] 0.016487137 0.032974273 0.9835129
[21,] 0.013184098 0.026368195 0.9868159
[22,] 0.010815425 0.021630849 0.9891846
[23,] 0.009652545 0.019305090 0.9903475
[24,] 0.008796657 0.017593313 0.9912033
[25,] 0.008025570 0.016051139 0.9919744
[26,] 0.009593746 0.019187492 0.9904063
[27,] 0.018655352 0.037310704 0.9813446
[28,] 0.018581176 0.037162351 0.9814188
[29,] 0.027400348 0.054800697 0.9725997
[30,] 0.031569883 0.063139767 0.9684301
[31,] 0.029323660 0.058647320 0.9706763
[32,] 0.031558177 0.063116355 0.9684418
[33,] 0.032773271 0.065546541 0.9672267
[34,] 0.033048343 0.066096686 0.9669517
[35,] 0.035113856 0.070227712 0.9648861
[36,] 0.051929616 0.103859233 0.9480704
[37,] 0.130536808 0.261073615 0.8694632
[38,] 0.245628895 0.491257789 0.7543711
[39,] 0.500283293 0.999433415 0.4997167
[40,] 0.593865368 0.812269264 0.4061346
[41,] 0.819559737 0.360880527 0.1804403
[42,] 0.802068837 0.395862326 0.1979312
[43,] 0.763425354 0.473149293 0.2365746
[44,] 0.701867273 0.596265453 0.2981327
[45,] 0.657373912 0.685252175 0.3426261
[46,] 0.635242818 0.729514365 0.3647572
[47,] 0.526628762 0.946742475 0.4733712
[48,] 0.461246036 0.922492072 0.5387540
[49,] 0.509130249 0.981739501 0.4908698
[50,] 0.682875794 0.634248413 0.3171242
[51,] 0.697700042 0.604599916 0.3023000
> postscript(file="/var/www/html/freestat/rcomp/tmp/17wms1229852811.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> plot(x[,1], type='l', main='Actuals and Interpolation', ylab='value of Actuals and Interpolation (dots)', xlab='time or index')
> points(x[,1]-mysum$resid)
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/2b69w1229852811.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> plot(mysum$resid, type='b', pch=19, main='Residuals', ylab='value of Residuals', xlab='time or index')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/3cfly1229852811.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> hist(mysum$resid, main='Residual Histogram', xlab='values of Residuals')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/4wo361229852811.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> densityplot(~mysum$resid,col='black',main='Residual Density Plot', xlab='values of Residuals')
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/5v24n1229852811.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 = 84
Frequency = 1
1 2 3 4 5 6
14699.3357 19224.3583 12617.1009 11533.8673 17091.4595 1974.6025
7 8 9 10 11 12
-7991.7467 8812.2889 -3814.0898 -5966.0620 -9485.3619 -2813.1406
13 14 15 16 17 18
-14429.9297 -10993.6057 -8436.5920 -9669.9895 -5664.1525 -8479.7188
19 20 21 22 23 24
-12579.2433 -11686.7755 -10787.6708 -5445.6363 -6996.7156 -3020.6839
25 26 27 28 29 30
-4632.7668 -8929.4692 -8715.0597 -4461.2420 -4175.4685 -8437.1608
31 32 33 34 35 36
-6221.8243 -185.5311 -6591.1242 -1586.9764 -762.5406 3761.6928
37 38 39 40 41 42
3324.1128 -858.6334 -11149.1309 -4145.4393 -4468.3577 -1311.1565
43 44 45 46 47 48
6622.3190 4673.0511 9951.6361 5803.4999 -442.1020 10463.2514
49 50 51 52 53 54
2816.6768 -832.3526 3132.0576 -2832.2205 509.9434 770.6476
55 56 57 58 59 60
9369.8972 13128.5221 14201.7544 19564.0835 23427.4959 17502.6772
61 62 63 64 65 66
15636.5821 11805.2951 16087.2461 9435.1817 3719.3571 14073.2893
67 68 69 70 71 72
20716.3124 8993.4017 11602.3198 14653.2411 13890.2629 11887.7720
73 74 75 76 77 78
-17414.0110 -9415.5926 -3535.6220 139.8423 -7012.7814 1409.4967
79 80 81 82 83 84
-9915.7143 -23734.9573 -14562.8255 -27022.1497 -19631.0387 -37781.5689
> postscript(file="/var/www/html/freestat/rcomp/tmp/6yhw01229852811.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 = 84
Frequency = 1
lag(myerror, k = 1) myerror
0 14699.3357 NA
1 19224.3583 14699.3357
2 12617.1009 19224.3583
3 11533.8673 12617.1009
4 17091.4595 11533.8673
5 1974.6025 17091.4595
6 -7991.7467 1974.6025
7 8812.2889 -7991.7467
8 -3814.0898 8812.2889
9 -5966.0620 -3814.0898
10 -9485.3619 -5966.0620
11 -2813.1406 -9485.3619
12 -14429.9297 -2813.1406
13 -10993.6057 -14429.9297
14 -8436.5920 -10993.6057
15 -9669.9895 -8436.5920
16 -5664.1525 -9669.9895
17 -8479.7188 -5664.1525
18 -12579.2433 -8479.7188
19 -11686.7755 -12579.2433
20 -10787.6708 -11686.7755
21 -5445.6363 -10787.6708
22 -6996.7156 -5445.6363
23 -3020.6839 -6996.7156
24 -4632.7668 -3020.6839
25 -8929.4692 -4632.7668
26 -8715.0597 -8929.4692
27 -4461.2420 -8715.0597
28 -4175.4685 -4461.2420
29 -8437.1608 -4175.4685
30 -6221.8243 -8437.1608
31 -185.5311 -6221.8243
32 -6591.1242 -185.5311
33 -1586.9764 -6591.1242
34 -762.5406 -1586.9764
35 3761.6928 -762.5406
36 3324.1128 3761.6928
37 -858.6334 3324.1128
38 -11149.1309 -858.6334
39 -4145.4393 -11149.1309
40 -4468.3577 -4145.4393
41 -1311.1565 -4468.3577
42 6622.3190 -1311.1565
43 4673.0511 6622.3190
44 9951.6361 4673.0511
45 5803.4999 9951.6361
46 -442.1020 5803.4999
47 10463.2514 -442.1020
48 2816.6768 10463.2514
49 -832.3526 2816.6768
50 3132.0576 -832.3526
51 -2832.2205 3132.0576
52 509.9434 -2832.2205
53 770.6476 509.9434
54 9369.8972 770.6476
55 13128.5221 9369.8972
56 14201.7544 13128.5221
57 19564.0835 14201.7544
58 23427.4959 19564.0835
59 17502.6772 23427.4959
60 15636.5821 17502.6772
61 11805.2951 15636.5821
62 16087.2461 11805.2951
63 9435.1817 16087.2461
64 3719.3571 9435.1817
65 14073.2893 3719.3571
66 20716.3124 14073.2893
67 8993.4017 20716.3124
68 11602.3198 8993.4017
69 14653.2411 11602.3198
70 13890.2629 14653.2411
71 11887.7720 13890.2629
72 -17414.0110 11887.7720
73 -9415.5926 -17414.0110
74 -3535.6220 -9415.5926
75 139.8423 -3535.6220
76 -7012.7814 139.8423
77 1409.4967 -7012.7814
78 -9915.7143 1409.4967
79 -23734.9573 -9915.7143
80 -14562.8255 -23734.9573
81 -27022.1497 -14562.8255
82 -19631.0387 -27022.1497
83 -37781.5689 -19631.0387
84 NA -37781.5689
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 19224.3583 14699.3357
[2,] 12617.1009 19224.3583
[3,] 11533.8673 12617.1009
[4,] 17091.4595 11533.8673
[5,] 1974.6025 17091.4595
[6,] -7991.7467 1974.6025
[7,] 8812.2889 -7991.7467
[8,] -3814.0898 8812.2889
[9,] -5966.0620 -3814.0898
[10,] -9485.3619 -5966.0620
[11,] -2813.1406 -9485.3619
[12,] -14429.9297 -2813.1406
[13,] -10993.6057 -14429.9297
[14,] -8436.5920 -10993.6057
[15,] -9669.9895 -8436.5920
[16,] -5664.1525 -9669.9895
[17,] -8479.7188 -5664.1525
[18,] -12579.2433 -8479.7188
[19,] -11686.7755 -12579.2433
[20,] -10787.6708 -11686.7755
[21,] -5445.6363 -10787.6708
[22,] -6996.7156 -5445.6363
[23,] -3020.6839 -6996.7156
[24,] -4632.7668 -3020.6839
[25,] -8929.4692 -4632.7668
[26,] -8715.0597 -8929.4692
[27,] -4461.2420 -8715.0597
[28,] -4175.4685 -4461.2420
[29,] -8437.1608 -4175.4685
[30,] -6221.8243 -8437.1608
[31,] -185.5311 -6221.8243
[32,] -6591.1242 -185.5311
[33,] -1586.9764 -6591.1242
[34,] -762.5406 -1586.9764
[35,] 3761.6928 -762.5406
[36,] 3324.1128 3761.6928
[37,] -858.6334 3324.1128
[38,] -11149.1309 -858.6334
[39,] -4145.4393 -11149.1309
[40,] -4468.3577 -4145.4393
[41,] -1311.1565 -4468.3577
[42,] 6622.3190 -1311.1565
[43,] 4673.0511 6622.3190
[44,] 9951.6361 4673.0511
[45,] 5803.4999 9951.6361
[46,] -442.1020 5803.4999
[47,] 10463.2514 -442.1020
[48,] 2816.6768 10463.2514
[49,] -832.3526 2816.6768
[50,] 3132.0576 -832.3526
[51,] -2832.2205 3132.0576
[52,] 509.9434 -2832.2205
[53,] 770.6476 509.9434
[54,] 9369.8972 770.6476
[55,] 13128.5221 9369.8972
[56,] 14201.7544 13128.5221
[57,] 19564.0835 14201.7544
[58,] 23427.4959 19564.0835
[59,] 17502.6772 23427.4959
[60,] 15636.5821 17502.6772
[61,] 11805.2951 15636.5821
[62,] 16087.2461 11805.2951
[63,] 9435.1817 16087.2461
[64,] 3719.3571 9435.1817
[65,] 14073.2893 3719.3571
[66,] 20716.3124 14073.2893
[67,] 8993.4017 20716.3124
[68,] 11602.3198 8993.4017
[69,] 14653.2411 11602.3198
[70,] 13890.2629 14653.2411
[71,] 11887.7720 13890.2629
[72,] -17414.0110 11887.7720
[73,] -9415.5926 -17414.0110
[74,] -3535.6220 -9415.5926
[75,] 139.8423 -3535.6220
[76,] -7012.7814 139.8423
[77,] 1409.4967 -7012.7814
[78,] -9915.7143 1409.4967
[79,] -23734.9573 -9915.7143
[80,] -14562.8255 -23734.9573
[81,] -27022.1497 -14562.8255
[82,] -19631.0387 -27022.1497
[83,] -37781.5689 -19631.0387
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 19224.3583 14699.3357
2 12617.1009 19224.3583
3 11533.8673 12617.1009
4 17091.4595 11533.8673
5 1974.6025 17091.4595
6 -7991.7467 1974.6025
7 8812.2889 -7991.7467
8 -3814.0898 8812.2889
9 -5966.0620 -3814.0898
10 -9485.3619 -5966.0620
11 -2813.1406 -9485.3619
12 -14429.9297 -2813.1406
13 -10993.6057 -14429.9297
14 -8436.5920 -10993.6057
15 -9669.9895 -8436.5920
16 -5664.1525 -9669.9895
17 -8479.7188 -5664.1525
18 -12579.2433 -8479.7188
19 -11686.7755 -12579.2433
20 -10787.6708 -11686.7755
21 -5445.6363 -10787.6708
22 -6996.7156 -5445.6363
23 -3020.6839 -6996.7156
24 -4632.7668 -3020.6839
25 -8929.4692 -4632.7668
26 -8715.0597 -8929.4692
27 -4461.2420 -8715.0597
28 -4175.4685 -4461.2420
29 -8437.1608 -4175.4685
30 -6221.8243 -8437.1608
31 -185.5311 -6221.8243
32 -6591.1242 -185.5311
33 -1586.9764 -6591.1242
34 -762.5406 -1586.9764
35 3761.6928 -762.5406
36 3324.1128 3761.6928
37 -858.6334 3324.1128
38 -11149.1309 -858.6334
39 -4145.4393 -11149.1309
40 -4468.3577 -4145.4393
41 -1311.1565 -4468.3577
42 6622.3190 -1311.1565
43 4673.0511 6622.3190
44 9951.6361 4673.0511
45 5803.4999 9951.6361
46 -442.1020 5803.4999
47 10463.2514 -442.1020
48 2816.6768 10463.2514
49 -832.3526 2816.6768
50 3132.0576 -832.3526
51 -2832.2205 3132.0576
52 509.9434 -2832.2205
53 770.6476 509.9434
54 9369.8972 770.6476
55 13128.5221 9369.8972
56 14201.7544 13128.5221
57 19564.0835 14201.7544
58 23427.4959 19564.0835
59 17502.6772 23427.4959
60 15636.5821 17502.6772
61 11805.2951 15636.5821
62 16087.2461 11805.2951
63 9435.1817 16087.2461
64 3719.3571 9435.1817
65 14073.2893 3719.3571
66 20716.3124 14073.2893
67 8993.4017 20716.3124
68 11602.3198 8993.4017
69 14653.2411 11602.3198
70 13890.2629 14653.2411
71 11887.7720 13890.2629
72 -17414.0110 11887.7720
73 -9415.5926 -17414.0110
74 -3535.6220 -9415.5926
75 139.8423 -3535.6220
76 -7012.7814 139.8423
77 1409.4967 -7012.7814
78 -9915.7143 1409.4967
79 -23734.9573 -9915.7143
80 -14562.8255 -23734.9573
81 -27022.1497 -14562.8255
82 -19631.0387 -27022.1497
83 -37781.5689 -19631.0387
> plot(z,main=paste('Residual Lag plot, lowess, and regression line'), ylab='values of Residuals', xlab='lagged values of Residuals')
> lines(lowess(z))
> abline(lm(z))
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/74tqz1229852811.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> acf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Autocorrelation Function')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/8v8m21229852811.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> pacf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Partial Autocorrelation Function')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/9huy81229852811.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> opar <- par(mfrow = c(2,2), oma = c(0, 0, 1.1, 0))
> plot(mylm, las = 1, sub='Residual Diagnostics')
> par(opar)
> dev.off()
null device
1
> if (n > n25) {
+ postscript(file="/var/www/html/freestat/rcomp/tmp/1001hm1229852811.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
+ plot(kp3:nmkm3,gqarr[,2], main='Goldfeld-Quandt test',ylab='2-sided p-value',xlab='breakpoint')
+ grid()
+ dev.off()
+ }
null device
1
>
> #Note: the /var/www/html/freestat/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/html/freestat/rcomp/createtable")
>
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Estimated Regression Equation', 1, TRUE)
> a<-table.row.end(a)
> myeq <- colnames(x)[1]
> myeq <- paste(myeq, '[t] = ', sep='')
> for (i in 1:k){
+ if (mysum$coefficients[i,1] > 0) myeq <- paste(myeq, '+', '')
+ myeq <- paste(myeq, mysum$coefficients[i,1], sep=' ')
+ if (rownames(mysum$coefficients)[i] != '(Intercept)') {
+ myeq <- paste(myeq, rownames(mysum$coefficients)[i], sep='')
+ if (rownames(mysum$coefficients)[i] != 't') myeq <- paste(myeq, '[t]', sep='')
+ }
+ }
> myeq <- paste(myeq, ' + e[t]')
> a<-table.row.start(a)
> a<-table.element(a, myeq)
> a<-table.row.end(a)
> a<-table.end(a)
> table.save(a,file="/var/www/html/freestat/rcomp/tmp/11eqoq1229852811.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a,hyperlink('http://www.xycoon.com/ols1.htm','Multiple Linear Regression - Ordinary Least Squares',''), 6, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a,'Variable',header=TRUE)
> a<-table.element(a,'Parameter',header=TRUE)
> a<-table.element(a,'S.D.',header=TRUE)
> a<-table.element(a,'T-STAT
H0: parameter = 0',header=TRUE)
> a<-table.element(a,'2-tail p-value',header=TRUE)
> a<-table.element(a,'1-tail p-value',header=TRUE)
> a<-table.row.end(a)
> for (i in 1:k){
+ a<-table.row.start(a)
+ a<-table.element(a,rownames(mysum$coefficients)[i],header=TRUE)
+ a<-table.element(a,mysum$coefficients[i,1])
+ a<-table.element(a, round(mysum$coefficients[i,2],6))
+ a<-table.element(a, round(mysum$coefficients[i,3],4))
+ a<-table.element(a, round(mysum$coefficients[i,4],6))
+ a<-table.element(a, round(mysum$coefficients[i,4]/2,6))
+ a<-table.row.end(a)
+ }
> a<-table.end(a)
> table.save(a,file="/var/www/html/freestat/rcomp/tmp/12rnm01229852811.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Regression Statistics', 2, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple R',1,TRUE)
> a<-table.element(a, sqrt(mysum$r.squared))
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'R-squared',1,TRUE)
> a<-table.element(a, mysum$r.squared)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Adjusted R-squared',1,TRUE)
> a<-table.element(a, mysum$adj.r.squared)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (value)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[1])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (DF numerator)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[2])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (DF denominator)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[3])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'p-value',1,TRUE)
> a<-table.element(a, 1-pf(mysum$fstatistic[1],mysum$fstatistic[2],mysum$fstatistic[3]))
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Residual Statistics', 2, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Residual Standard Deviation',1,TRUE)
> a<-table.element(a, mysum$sigma)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Sum Squared Residuals',1,TRUE)
> a<-table.element(a, sum(myerror*myerror))
> a<-table.row.end(a)
> a<-table.end(a)
> table.save(a,file="/var/www/html/freestat/rcomp/tmp/137vil1229852811.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Actuals, Interpolation, and Residuals', 4, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Time or Index', 1, TRUE)
> a<-table.element(a, 'Actuals', 1, TRUE)
> a<-table.element(a, 'Interpolation
Forecast', 1, TRUE)
> a<-table.element(a, 'Residuals
Prediction Error', 1, TRUE)
> a<-table.row.end(a)
> for (i in 1:n) {
+ a<-table.row.start(a)
+ a<-table.element(a,i, 1, TRUE)
+ a<-table.element(a,x[i])
+ a<-table.element(a,x[i]-mysum$resid[i])
+ a<-table.element(a,mysum$resid[i])
+ a<-table.row.end(a)
+ }
> a<-table.end(a)
> table.save(a,file="/var/www/html/freestat/rcomp/tmp/14vafp1229852811.tab")
> if (n > n25) {
+ a<-table.start()
+ a<-table.row.start(a)
+ a<-table.element(a,'Goldfeld-Quandt test for Heteroskedasticity',4,TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'p-values',header=TRUE)
+ a<-table.element(a,'Alternative Hypothesis',3,header=TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'breakpoint index',header=TRUE)
+ a<-table.element(a,'greater',header=TRUE)
+ a<-table.element(a,'2-sided',header=TRUE)
+ a<-table.element(a,'less',header=TRUE)
+ a<-table.row.end(a)
+ for (mypoint in kp3:nmkm3) {
+ a<-table.row.start(a)
+ a<-table.element(a,mypoint,header=TRUE)
+ a<-table.element(a,gqarr[mypoint-kp3+1,1])
+ a<-table.element(a,gqarr[mypoint-kp3+1,2])
+ a<-table.element(a,gqarr[mypoint-kp3+1,3])
+ a<-table.row.end(a)
+ }
+ a<-table.end(a)
+ table.save(a,file="/var/www/html/freestat/rcomp/tmp/1531av1229852811.tab")
+ a<-table.start()
+ a<-table.row.start(a)
+ a<-table.element(a,'Meta Analysis of Goldfeld-Quandt test for Heteroskedasticity',4,TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'Description',header=TRUE)
+ a<-table.element(a,'# significant tests',header=TRUE)
+ a<-table.element(a,'% significant tests',header=TRUE)
+ a<-table.element(a,'OK/NOK',header=TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'1% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant1)
+ a<-table.element(a,numsignificant1/numgqtests)
+ if (numsignificant1/numgqtests < 0.01) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'5% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant5)
+ a<-table.element(a,numsignificant5/numgqtests)
+ if (numsignificant5/numgqtests < 0.05) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'10% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant10)
+ a<-table.element(a,numsignificant10/numgqtests)
+ if (numsignificant10/numgqtests < 0.1) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.end(a)
+ table.save(a,file="/var/www/html/freestat/rcomp/tmp/16ebvw1229852811.tab")
+ }
>
> system("convert tmp/17wms1229852811.ps tmp/17wms1229852811.png")
> system("convert tmp/2b69w1229852811.ps tmp/2b69w1229852811.png")
> system("convert tmp/3cfly1229852811.ps tmp/3cfly1229852811.png")
> system("convert tmp/4wo361229852811.ps tmp/4wo361229852811.png")
> system("convert tmp/5v24n1229852811.ps tmp/5v24n1229852811.png")
> system("convert tmp/6yhw01229852811.ps tmp/6yhw01229852811.png")
> system("convert tmp/74tqz1229852811.ps tmp/74tqz1229852811.png")
> system("convert tmp/8v8m21229852811.ps tmp/8v8m21229852811.png")
> system("convert tmp/9huy81229852811.ps tmp/9huy81229852811.png")
> system("convert tmp/1001hm1229852811.ps tmp/1001hm1229852811.png")
>
>
> proc.time()
user system elapsed
4.083 2.567 4.658