R version 2.6.1 (2007-11-26)
Copyright (C) 2007 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(13328,12873,14000,13477,14237,13674,13529,14058,12975,14326,14008,16193,14483,14011,15057,14884,15414,14440,14900,15074,14442,15307,14938,17193,15528,14765,15838,15723,16150,15486,15986,15983,15692,16490,15686,18897,16316,15636,17163,16534,16518,16375,16290,16352,15943,16362,16393,19051,16747,16320,17910,16961,17480,17049,16879,17473,16998,17307,17418,20169,17871,17226,19062,17804,19100,18522,18060,18869,18127,18871,18890,21263,19547,18450,20254,19240,20216,19420,19415,20018,18652,19978,19509,21971),dim=c(1,84),dimnames=list(c('HPC'),1:84))
> y <- array(NA,dim=c(1,84),dimnames=list(c('HPC'),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
> 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
HPC M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 13328 1 0 0 0 0 0 0 0 0 0 0 1
2 12873 0 1 0 0 0 0 0 0 0 0 0 2
3 14000 0 0 1 0 0 0 0 0 0 0 0 3
4 13477 0 0 0 1 0 0 0 0 0 0 0 4
5 14237 0 0 0 0 1 0 0 0 0 0 0 5
6 13674 0 0 0 0 0 1 0 0 0 0 0 6
7 13529 0 0 0 0 0 0 1 0 0 0 0 7
8 14058 0 0 0 0 0 0 0 1 0 0 0 8
9 12975 0 0 0 0 0 0 0 0 1 0 0 9
10 14326 0 0 0 0 0 0 0 0 0 1 0 10
11 14008 0 0 0 0 0 0 0 0 0 0 1 11
12 16193 0 0 0 0 0 0 0 0 0 0 0 12
13 14483 1 0 0 0 0 0 0 0 0 0 0 13
14 14011 0 1 0 0 0 0 0 0 0 0 0 14
15 15057 0 0 1 0 0 0 0 0 0 0 0 15
16 14884 0 0 0 1 0 0 0 0 0 0 0 16
17 15414 0 0 0 0 1 0 0 0 0 0 0 17
18 14440 0 0 0 0 0 1 0 0 0 0 0 18
19 14900 0 0 0 0 0 0 1 0 0 0 0 19
20 15074 0 0 0 0 0 0 0 1 0 0 0 20
21 14442 0 0 0 0 0 0 0 0 1 0 0 21
22 15307 0 0 0 0 0 0 0 0 0 1 0 22
23 14938 0 0 0 0 0 0 0 0 0 0 1 23
24 17193 0 0 0 0 0 0 0 0 0 0 0 24
25 15528 1 0 0 0 0 0 0 0 0 0 0 25
26 14765 0 1 0 0 0 0 0 0 0 0 0 26
27 15838 0 0 1 0 0 0 0 0 0 0 0 27
28 15723 0 0 0 1 0 0 0 0 0 0 0 28
29 16150 0 0 0 0 1 0 0 0 0 0 0 29
30 15486 0 0 0 0 0 1 0 0 0 0 0 30
31 15986 0 0 0 0 0 0 1 0 0 0 0 31
32 15983 0 0 0 0 0 0 0 1 0 0 0 32
33 15692 0 0 0 0 0 0 0 0 1 0 0 33
34 16490 0 0 0 0 0 0 0 0 0 1 0 34
35 15686 0 0 0 0 0 0 0 0 0 0 1 35
36 18897 0 0 0 0 0 0 0 0 0 0 0 36
37 16316 1 0 0 0 0 0 0 0 0 0 0 37
38 15636 0 1 0 0 0 0 0 0 0 0 0 38
39 17163 0 0 1 0 0 0 0 0 0 0 0 39
40 16534 0 0 0 1 0 0 0 0 0 0 0 40
41 16518 0 0 0 0 1 0 0 0 0 0 0 41
42 16375 0 0 0 0 0 1 0 0 0 0 0 42
43 16290 0 0 0 0 0 0 1 0 0 0 0 43
44 16352 0 0 0 0 0 0 0 1 0 0 0 44
45 15943 0 0 0 0 0 0 0 0 1 0 0 45
46 16362 0 0 0 0 0 0 0 0 0 1 0 46
47 16393 0 0 0 0 0 0 0 0 0 0 1 47
48 19051 0 0 0 0 0 0 0 0 0 0 0 48
49 16747 1 0 0 0 0 0 0 0 0 0 0 49
50 16320 0 1 0 0 0 0 0 0 0 0 0 50
51 17910 0 0 1 0 0 0 0 0 0 0 0 51
52 16961 0 0 0 1 0 0 0 0 0 0 0 52
53 17480 0 0 0 0 1 0 0 0 0 0 0 53
54 17049 0 0 0 0 0 1 0 0 0 0 0 54
55 16879 0 0 0 0 0 0 1 0 0 0 0 55
56 17473 0 0 0 0 0 0 0 1 0 0 0 56
57 16998 0 0 0 0 0 0 0 0 1 0 0 57
58 17307 0 0 0 0 0 0 0 0 0 1 0 58
59 17418 0 0 0 0 0 0 0 0 0 0 1 59
60 20169 0 0 0 0 0 0 0 0 0 0 0 60
61 17871 1 0 0 0 0 0 0 0 0 0 0 61
62 17226 0 1 0 0 0 0 0 0 0 0 0 62
63 19062 0 0 1 0 0 0 0 0 0 0 0 63
64 17804 0 0 0 1 0 0 0 0 0 0 0 64
65 19100 0 0 0 0 1 0 0 0 0 0 0 65
66 18522 0 0 0 0 0 1 0 0 0 0 0 66
67 18060 0 0 0 0 0 0 1 0 0 0 0 67
68 18869 0 0 0 0 0 0 0 1 0 0 0 68
69 18127 0 0 0 0 0 0 0 0 1 0 0 69
70 18871 0 0 0 0 0 0 0 0 0 1 0 70
71 18890 0 0 0 0 0 0 0 0 0 0 1 71
72 21263 0 0 0 0 0 0 0 0 0 0 0 72
73 19547 1 0 0 0 0 0 0 0 0 0 0 73
74 18450 0 1 0 0 0 0 0 0 0 0 0 74
75 20254 0 0 1 0 0 0 0 0 0 0 0 75
76 19240 0 0 0 1 0 0 0 0 0 0 0 76
77 20216 0 0 0 0 1 0 0 0 0 0 0 77
78 19420 0 0 0 0 0 1 0 0 0 0 0 78
79 19415 0 0 0 0 0 0 1 0 0 0 0 79
80 20018 0 0 0 0 0 0 0 1 0 0 0 80
81 18652 0 0 0 0 0 0 0 0 1 0 0 81
82 19978 0 0 0 0 0 0 0 0 0 1 0 82
83 19509 0 0 0 0 0 0 0 0 0 0 1 83
84 21971 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) M1 M2 M3 M4 M5
15515.92 -2132.84 -2859.02 -1507.78 -2251.39 -1687.43
M6 M7 M8 M9 M10 M11
-2357.90 -2422.37 -2104.70 -2896.45 -2143.92 -2478.67
t
77.75
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-586.7 -211.2 32.0 181.5 581.9
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 15515.917 129.278 120.020 < 2e-16 ***
M1 -2132.841 159.024 -13.412 < 2e-16 ***
M2 -2859.024 158.904 -17.992 < 2e-16 ***
M3 -1507.779 158.796 -9.495 2.88e-14 ***
M4 -2251.391 158.699 -14.187 < 2e-16 ***
M5 -1687.431 158.613 -10.639 2.40e-16 ***
M6 -2357.900 158.538 -14.873 < 2e-16 ***
M7 -2422.369 158.476 -15.285 < 2e-16 ***
M8 -2104.695 158.424 -13.285 < 2e-16 ***
M9 -2896.450 158.384 -18.288 < 2e-16 ***
M10 -2143.919 158.355 -13.539 < 2e-16 ***
M11 -2478.674 158.338 -15.654 < 2e-16 ***
t 77.755 1.347 57.740 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 296.2 on 71 degrees of freedom
Multiple R-Squared: 0.9827, Adjusted R-squared: 0.9798
F-statistic: 335.8 on 12 and 71 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.0480796324 0.0961592648 0.9519204
[2,] 0.0132432467 0.0264864934 0.9867568
[3,] 0.0432107824 0.0864215648 0.9567892
[4,] 0.0315084969 0.0630169938 0.9684915
[5,] 0.0152933905 0.0305867811 0.9847066
[6,] 0.0158812132 0.0317624263 0.9841188
[7,] 0.0096668261 0.0193336523 0.9903332
[8,] 0.0063019844 0.0126039687 0.9936980
[9,] 0.0031909454 0.0063818908 0.9968091
[10,] 0.0014335262 0.0028670524 0.9985665
[11,] 0.0021211010 0.0042422019 0.9978789
[12,] 0.0023861570 0.0047723139 0.9976138
[13,] 0.0013142002 0.0026284005 0.9986858
[14,] 0.0010446325 0.0020892650 0.9989554
[15,] 0.0005141695 0.0010283389 0.9994858
[16,] 0.0007761228 0.0015522455 0.9992239
[17,] 0.0004520142 0.0009040284 0.9995480
[18,] 0.0023910490 0.0047820980 0.9976090
[19,] 0.0048158827 0.0096317654 0.9951841
[20,] 0.0068661724 0.0137323448 0.9931338
[21,] 0.1022724671 0.2045449343 0.8977275
[22,] 0.1133385592 0.2266771184 0.8866614
[23,] 0.1597194998 0.3194389996 0.8402805
[24,] 0.1391865741 0.2783731481 0.8608134
[25,] 0.2915853738 0.5831707476 0.7084146
[26,] 0.5844744742 0.8310510517 0.4155255
[27,] 0.5646397968 0.8707204064 0.4353602
[28,] 0.7026217470 0.5947565060 0.2973783
[29,] 0.7823532289 0.4352935423 0.2176468
[30,] 0.7976629207 0.4046741586 0.2023371
[31,] 0.8983544045 0.2032911910 0.1016456
[32,] 0.8768903564 0.2462192872 0.1231096
[33,] 0.8586538490 0.2826923020 0.1413462
[34,] 0.8571553959 0.2856892081 0.1428446
[35,] 0.8348847469 0.3302305063 0.1651153
[36,] 0.7856957300 0.4286085400 0.2143043
[37,] 0.7729548268 0.4540903464 0.2270452
[38,] 0.8095783383 0.3808433235 0.1904217
[39,] 0.7708311471 0.4583377059 0.2291689
[40,] 0.7537279227 0.4925441545 0.2462721
[41,] 0.7303457188 0.5393085624 0.2696543
[42,] 0.7019123474 0.5961753052 0.2980877
[43,] 0.7867422343 0.4265155313 0.2132578
[44,] 0.7314442686 0.5371114629 0.2685557
[45,] 0.6697964077 0.6604071846 0.3302036
[46,] 0.8007811678 0.3984376644 0.1992188
[47,] 0.7441795762 0.5116408477 0.2558204
[48,] 0.7041095070 0.5917809861 0.2958905
[49,] 0.7764032321 0.4471935358 0.2235968
[50,] 0.7306154255 0.5387691490 0.2693846
[51,] 0.6287885172 0.7424229657 0.3712115
[52,] 0.7246310493 0.5507379013 0.2753690
[53,] 0.7528527866 0.4942944268 0.2471472
> postscript(file="/var/www/html/rcomp/tmp/1oiwj1204983648.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/2c4xh1204983648.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/3b8b61204983649.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/4l4df1204983649.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/52n8o1204983649.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
-132.830357 60.598214 -241.401786 -98.544643 19.741071 49.455357
7 8 9 10 11 12
-108.830357 24.741071 -344.258929 176.455357 115.455357 -255.973214
13 14 15 16 17 18
89.113095 265.541667 -117.458333 375.398810 263.684524 -117.601190
19 20 21 22 23 24
329.113095 107.684524 189.684524 224.398810 112.398810 -189.029762
25 26 27 28 29 30
201.056548 86.485119 -269.514881 281.342262 66.627976 -4.657738
31 32 33 34 35 36
482.056548 83.627976 506.627976 474.342262 -72.657738 581.913690
37 38 39 40 41 42
56.000000 24.428571 122.428571 159.285714 -498.428571 -48.714286
43 44 45 46 47 48
-147.000000 -480.428571 -175.428571 -586.714286 -298.714286 -197.142857
49 50 51 52 53 54
-446.056548 -224.627976 -63.627976 -346.770833 -469.485119 -307.770833
55 56 57 58 59 60
-491.056548 -292.485119 -53.485119 -574.770833 -206.770833 -12.199405
61 62 63 64 65 66
-255.113095 -251.684524 155.315476 -436.827381 217.458333 232.172619
67 68 69 70 71 72
-243.113095 170.458333 142.458333 56.172619 332.172619 148.744048
73 74 75 76 77 78
487.830357 39.258929 414.258929 66.116071 400.401786 197.116071
79 80 81 82 83 84
178.830357 386.401786 -265.598214 230.116071 18.116071 -76.312500
> postscript(file="/var/www/html/rcomp/tmp/6io7c1204983649.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 -132.830357 NA
1 60.598214 -132.830357
2 -241.401786 60.598214
3 -98.544643 -241.401786
4 19.741071 -98.544643
5 49.455357 19.741071
6 -108.830357 49.455357
7 24.741071 -108.830357
8 -344.258929 24.741071
9 176.455357 -344.258929
10 115.455357 176.455357
11 -255.973214 115.455357
12 89.113095 -255.973214
13 265.541667 89.113095
14 -117.458333 265.541667
15 375.398810 -117.458333
16 263.684524 375.398810
17 -117.601190 263.684524
18 329.113095 -117.601190
19 107.684524 329.113095
20 189.684524 107.684524
21 224.398810 189.684524
22 112.398810 224.398810
23 -189.029762 112.398810
24 201.056548 -189.029762
25 86.485119 201.056548
26 -269.514881 86.485119
27 281.342262 -269.514881
28 66.627976 281.342262
29 -4.657738 66.627976
30 482.056548 -4.657738
31 83.627976 482.056548
32 506.627976 83.627976
33 474.342262 506.627976
34 -72.657738 474.342262
35 581.913690 -72.657738
36 56.000000 581.913690
37 24.428571 56.000000
38 122.428571 24.428571
39 159.285714 122.428571
40 -498.428571 159.285714
41 -48.714286 -498.428571
42 -147.000000 -48.714286
43 -480.428571 -147.000000
44 -175.428571 -480.428571
45 -586.714286 -175.428571
46 -298.714286 -586.714286
47 -197.142857 -298.714286
48 -446.056548 -197.142857
49 -224.627976 -446.056548
50 -63.627976 -224.627976
51 -346.770833 -63.627976
52 -469.485119 -346.770833
53 -307.770833 -469.485119
54 -491.056548 -307.770833
55 -292.485119 -491.056548
56 -53.485119 -292.485119
57 -574.770833 -53.485119
58 -206.770833 -574.770833
59 -12.199405 -206.770833
60 -255.113095 -12.199405
61 -251.684524 -255.113095
62 155.315476 -251.684524
63 -436.827381 155.315476
64 217.458333 -436.827381
65 232.172619 217.458333
66 -243.113095 232.172619
67 170.458333 -243.113095
68 142.458333 170.458333
69 56.172619 142.458333
70 332.172619 56.172619
71 148.744048 332.172619
72 487.830357 148.744048
73 39.258929 487.830357
74 414.258929 39.258929
75 66.116071 414.258929
76 400.401786 66.116071
77 197.116071 400.401786
78 178.830357 197.116071
79 386.401786 178.830357
80 -265.598214 386.401786
81 230.116071 -265.598214
82 18.116071 230.116071
83 -76.312500 18.116071
84 NA -76.312500
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 60.598214 -132.830357
[2,] -241.401786 60.598214
[3,] -98.544643 -241.401786
[4,] 19.741071 -98.544643
[5,] 49.455357 19.741071
[6,] -108.830357 49.455357
[7,] 24.741071 -108.830357
[8,] -344.258929 24.741071
[9,] 176.455357 -344.258929
[10,] 115.455357 176.455357
[11,] -255.973214 115.455357
[12,] 89.113095 -255.973214
[13,] 265.541667 89.113095
[14,] -117.458333 265.541667
[15,] 375.398810 -117.458333
[16,] 263.684524 375.398810
[17,] -117.601190 263.684524
[18,] 329.113095 -117.601190
[19,] 107.684524 329.113095
[20,] 189.684524 107.684524
[21,] 224.398810 189.684524
[22,] 112.398810 224.398810
[23,] -189.029762 112.398810
[24,] 201.056548 -189.029762
[25,] 86.485119 201.056548
[26,] -269.514881 86.485119
[27,] 281.342262 -269.514881
[28,] 66.627976 281.342262
[29,] -4.657738 66.627976
[30,] 482.056548 -4.657738
[31,] 83.627976 482.056548
[32,] 506.627976 83.627976
[33,] 474.342262 506.627976
[34,] -72.657738 474.342262
[35,] 581.913690 -72.657738
[36,] 56.000000 581.913690
[37,] 24.428571 56.000000
[38,] 122.428571 24.428571
[39,] 159.285714 122.428571
[40,] -498.428571 159.285714
[41,] -48.714286 -498.428571
[42,] -147.000000 -48.714286
[43,] -480.428571 -147.000000
[44,] -175.428571 -480.428571
[45,] -586.714286 -175.428571
[46,] -298.714286 -586.714286
[47,] -197.142857 -298.714286
[48,] -446.056548 -197.142857
[49,] -224.627976 -446.056548
[50,] -63.627976 -224.627976
[51,] -346.770833 -63.627976
[52,] -469.485119 -346.770833
[53,] -307.770833 -469.485119
[54,] -491.056548 -307.770833
[55,] -292.485119 -491.056548
[56,] -53.485119 -292.485119
[57,] -574.770833 -53.485119
[58,] -206.770833 -574.770833
[59,] -12.199405 -206.770833
[60,] -255.113095 -12.199405
[61,] -251.684524 -255.113095
[62,] 155.315476 -251.684524
[63,] -436.827381 155.315476
[64,] 217.458333 -436.827381
[65,] 232.172619 217.458333
[66,] -243.113095 232.172619
[67,] 170.458333 -243.113095
[68,] 142.458333 170.458333
[69,] 56.172619 142.458333
[70,] 332.172619 56.172619
[71,] 148.744048 332.172619
[72,] 487.830357 148.744048
[73,] 39.258929 487.830357
[74,] 414.258929 39.258929
[75,] 66.116071 414.258929
[76,] 400.401786 66.116071
[77,] 197.116071 400.401786
[78,] 178.830357 197.116071
[79,] 386.401786 178.830357
[80,] -265.598214 386.401786
[81,] 230.116071 -265.598214
[82,] 18.116071 230.116071
[83,] -76.312500 18.116071
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 60.598214 -132.830357
2 -241.401786 60.598214
3 -98.544643 -241.401786
4 19.741071 -98.544643
5 49.455357 19.741071
6 -108.830357 49.455357
7 24.741071 -108.830357
8 -344.258929 24.741071
9 176.455357 -344.258929
10 115.455357 176.455357
11 -255.973214 115.455357
12 89.113095 -255.973214
13 265.541667 89.113095
14 -117.458333 265.541667
15 375.398810 -117.458333
16 263.684524 375.398810
17 -117.601190 263.684524
18 329.113095 -117.601190
19 107.684524 329.113095
20 189.684524 107.684524
21 224.398810 189.684524
22 112.398810 224.398810
23 -189.029762 112.398810
24 201.056548 -189.029762
25 86.485119 201.056548
26 -269.514881 86.485119
27 281.342262 -269.514881
28 66.627976 281.342262
29 -4.657738 66.627976
30 482.056548 -4.657738
31 83.627976 482.056548
32 506.627976 83.627976
33 474.342262 506.627976
34 -72.657738 474.342262
35 581.913690 -72.657738
36 56.000000 581.913690
37 24.428571 56.000000
38 122.428571 24.428571
39 159.285714 122.428571
40 -498.428571 159.285714
41 -48.714286 -498.428571
42 -147.000000 -48.714286
43 -480.428571 -147.000000
44 -175.428571 -480.428571
45 -586.714286 -175.428571
46 -298.714286 -586.714286
47 -197.142857 -298.714286
48 -446.056548 -197.142857
49 -224.627976 -446.056548
50 -63.627976 -224.627976
51 -346.770833 -63.627976
52 -469.485119 -346.770833
53 -307.770833 -469.485119
54 -491.056548 -307.770833
55 -292.485119 -491.056548
56 -53.485119 -292.485119
57 -574.770833 -53.485119
58 -206.770833 -574.770833
59 -12.199405 -206.770833
60 -255.113095 -12.199405
61 -251.684524 -255.113095
62 155.315476 -251.684524
63 -436.827381 155.315476
64 217.458333 -436.827381
65 232.172619 217.458333
66 -243.113095 232.172619
67 170.458333 -243.113095
68 142.458333 170.458333
69 56.172619 142.458333
70 332.172619 56.172619
71 148.744048 332.172619
72 487.830357 148.744048
73 39.258929 487.830357
74 414.258929 39.258929
75 66.116071 414.258929
76 400.401786 66.116071
77 197.116071 400.401786
78 178.830357 197.116071
79 386.401786 178.830357
80 -265.598214 386.401786
81 230.116071 -265.598214
82 18.116071 230.116071
83 -76.312500 18.116071
> 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/7hsik1204983649.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/851161204983649.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/9k0e51204983649.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/10mr2f1204983649.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
> 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/11781m1204983649.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/121pcf1204983649.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/137jxr1204983649.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/14dudz1204983649.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/15sp9w1204983649.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/16ig3x1204983649.tab")
+ }
>
> system("convert tmp/1oiwj1204983648.ps tmp/1oiwj1204983648.png")
> system("convert tmp/2c4xh1204983648.ps tmp/2c4xh1204983648.png")
> system("convert tmp/3b8b61204983649.ps tmp/3b8b61204983649.png")
> system("convert tmp/4l4df1204983649.ps tmp/4l4df1204983649.png")
> system("convert tmp/52n8o1204983649.ps tmp/52n8o1204983649.png")
> system("convert tmp/6io7c1204983649.ps tmp/6io7c1204983649.png")
> system("convert tmp/7hsik1204983649.ps tmp/7hsik1204983649.png")
> system("convert tmp/851161204983649.ps tmp/851161204983649.png")
> system("convert tmp/9k0e51204983649.ps tmp/9k0e51204983649.png")
> system("convert tmp/10mr2f1204983649.ps tmp/10mr2f1204983649.png")
>
>
> proc.time()
user system elapsed
2.984 1.590 3.720