R version 2.9.0 (2009-04-17)
Copyright (C) 2009 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(9051,0,8823,0,8776,0,8255,0,7969,0,8758,0,8693,0,8271,0,7790,0,7769,0,8170,0,8209,0,9395,0,9260,0,9018,0,8501,0,8500,0,9649,0,9319,0,8830,0,8436,0,8169,0,8269,0,7945,0,9144,0,8770,0,8834,0,7837,0,7792,0,8616,0,8518,0,7940,0,7545,0,7531,0,7665,0,7599,0,8444,0,8549,0,7986,0,7335,0,7287,0,7870,0,7839,0,7327,0,7259,0,6964,0,7271,0,6956,0,7608,0,7692,0,7255,0,6804,0,6655,0,7341,0,7602,0,7086,0,6625,0,6272,0,6576,0,6491,0,7649,0,7400,0,6913,0,6532,0,6486,0,7295,0,7556,0,7088,1,6952,1,6773,1,6917,1,7371,1,8221,1,7953,1,8027,1,7287,1,8076,1,8933,1,9433,1,9479,1,9199,1,9469,1,10015,1,10999,1,13009,1,13699,1,13895,1,13248,1,13973,1,15095,1,15201,1,14823,1,14538,1,14547,1,14407,1),dim=c(2,95),dimnames=list(c('Y','X'),1:95))
> y <- array(NA,dim=c(2,95),dimnames=list(c('Y','X'),1:95))
> 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 = 'Do not include Seasonal Dummies'
> par1 = '1'
> #'GNU S' R Code compiled by R2WASP v. 1.0.44 ()
> #Author: Prof. Dr. P. Wessa
> #To cite this work: AUTHOR(S), (YEAR), YOUR SOFTWARE TITLE (vNUMBER) in Free Statistics Software (v$_version), Office for Research Development and Education, URL http://www.wessa.net/rwasp_YOURPAGE.wasp/
> #Source of accompanying publication: Office for Research, Development, and Education
> #Technical description: Write here your technical program description (don't use hard returns!)
> library(lattice)
> library(lmtest)
Loading required package: zoo
Attaching package: 'zoo'
The following object(s) are masked from package:base :
as.Date.numeric
> n25 <- 25 #minimum number of obs. for Goldfeld-Quandt test
> par1 <- as.numeric(par1)
> x <- t(y)
> k <- length(x[1,])
> n <- length(x[,1])
> x1 <- cbind(x[,par1], x[,1:k!=par1])
> mycolnames <- c(colnames(x)[par1], colnames(x)[1:k!=par1])
> colnames(x1) <- mycolnames #colnames(x)[par1]
> x <- x1
> if (par3 == 'First Differences'){
+ x2 <- array(0, dim=c(n-1,k), dimnames=list(1:(n-1), paste('(1-B)',colnames(x),sep='')))
+ for (i in 1:n-1) {
+ for (j in 1:k) {
+ x2[i,j] <- x[i+1,j] - x[i,j]
+ }
+ }
+ x <- x2
+ }
> if (par2 == 'Include Monthly Dummies'){
+ x2 <- array(0, dim=c(n,11), dimnames=list(1:n, paste('M', seq(1:11), sep ='')))
+ for (i in 1:11){
+ x2[seq(i,n,12),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> if (par2 == 'Include Quarterly Dummies'){
+ x2 <- array(0, dim=c(n,3), dimnames=list(1:n, paste('Q', seq(1:3), sep ='')))
+ for (i in 1:3){
+ x2[seq(i,n,4),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> k <- length(x[1,])
> if (par3 == 'Linear Trend'){
+ x <- cbind(x, c(1:n))
+ colnames(x)[k+1] <- 't'
+ }
> x
Y X t
1 9051 0 1
2 8823 0 2
3 8776 0 3
4 8255 0 4
5 7969 0 5
6 8758 0 6
7 8693 0 7
8 8271 0 8
9 7790 0 9
10 7769 0 10
11 8170 0 11
12 8209 0 12
13 9395 0 13
14 9260 0 14
15 9018 0 15
16 8501 0 16
17 8500 0 17
18 9649 0 18
19 9319 0 19
20 8830 0 20
21 8436 0 21
22 8169 0 22
23 8269 0 23
24 7945 0 24
25 9144 0 25
26 8770 0 26
27 8834 0 27
28 7837 0 28
29 7792 0 29
30 8616 0 30
31 8518 0 31
32 7940 0 32
33 7545 0 33
34 7531 0 34
35 7665 0 35
36 7599 0 36
37 8444 0 37
38 8549 0 38
39 7986 0 39
40 7335 0 40
41 7287 0 41
42 7870 0 42
43 7839 0 43
44 7327 0 44
45 7259 0 45
46 6964 0 46
47 7271 0 47
48 6956 0 48
49 7608 0 49
50 7692 0 50
51 7255 0 51
52 6804 0 52
53 6655 0 53
54 7341 0 54
55 7602 0 55
56 7086 0 56
57 6625 0 57
58 6272 0 58
59 6576 0 59
60 6491 0 60
61 7649 0 61
62 7400 0 62
63 6913 0 63
64 6532 0 64
65 6486 0 65
66 7295 0 66
67 7556 0 67
68 7088 1 68
69 6952 1 69
70 6773 1 70
71 6917 1 71
72 7371 1 72
73 8221 1 73
74 7953 1 74
75 8027 1 75
76 7287 1 76
77 8076 1 77
78 8933 1 78
79 9433 1 79
80 9479 1 80
81 9199 1 81
82 9469 1 82
83 10015 1 83
84 10999 1 84
85 13009 1 85
86 13699 1 86
87 13895 1 87
88 13248 1 88
89 13973 1 89
90 15095 1 90
91 15201 1 91
92 14823 1 92
93 14538 1 93
94 14547 1 94
95 14407 1 95
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X t
8100.425 3072.410 -6.228
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-3963.9 -847.7 -100.3 694.6 4594.9
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 8100.425 436.359 18.564 < 2e-16 ***
X 3072.410 665.087 4.620 1.25e-05 ***
t -6.228 11.058 -0.563 0.575
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1813 on 92 degrees of freedom
Multiple R-squared: 0.3364, Adjusted R-squared: 0.3219
F-statistic: 23.31 on 2 and 92 DF, p-value: 6.44e-09
> 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,] 1.524008e-02 3.048017e-02 0.984759917
[2,] 4.591339e-03 9.182679e-03 0.995408661
[3,] 8.610161e-04 1.722032e-03 0.999138984
[4,] 2.342349e-04 4.684697e-04 0.999765765
[5,] 4.410009e-05 8.820018e-05 0.999955900
[6,] 1.141872e-05 2.283743e-05 0.999988581
[7,] 3.092489e-06 6.184979e-06 0.999996908
[8,] 5.667609e-05 1.133522e-04 0.999943324
[9,] 5.237882e-05 1.047576e-04 0.999947621
[10,] 2.037361e-05 4.074722e-05 0.999979626
[11,] 5.757336e-06 1.151467e-05 0.999994243
[12,] 1.568554e-06 3.137108e-06 0.999998431
[13,] 1.935421e-06 3.870841e-06 0.999998065
[14,] 8.028075e-07 1.605615e-06 0.999999197
[15,] 2.448639e-07 4.897278e-07 0.999999755
[16,] 1.016087e-07 2.032174e-07 0.999999898
[17,] 5.764386e-08 1.152877e-07 0.999999942
[18,] 2.384932e-08 4.769863e-08 0.999999976
[19,] 1.465981e-08 2.931961e-08 0.999999985
[20,] 7.718554e-09 1.543711e-08 0.999999992
[21,] 2.683294e-09 5.366588e-09 0.999999997
[22,] 9.964548e-10 1.992910e-09 0.999999999
[23,] 9.313524e-10 1.862705e-09 0.999999999
[24,] 7.243291e-10 1.448658e-09 0.999999999
[25,] 3.010924e-10 6.021848e-10 1.000000000
[26,] 1.275049e-10 2.550098e-10 1.000000000
[27,] 7.606921e-11 1.521384e-10 1.000000000
[28,] 8.082282e-11 1.616456e-10 1.000000000
[29,] 6.927973e-11 1.385595e-10 1.000000000
[30,] 4.254292e-11 8.508584e-11 1.000000000
[31,] 2.650418e-11 5.300837e-11 1.000000000
[32,] 2.177318e-11 4.354637e-11 1.000000000
[33,] 2.529891e-11 5.059781e-11 1.000000000
[34,] 1.937655e-11 3.875310e-11 1.000000000
[35,] 2.420466e-11 4.840931e-11 1.000000000
[36,] 2.918292e-11 5.836583e-11 1.000000000
[37,] 3.081613e-11 6.163227e-11 1.000000000
[38,] 3.936135e-11 7.872271e-11 1.000000000
[39,] 5.741775e-11 1.148355e-10 1.000000000
[40,] 9.058980e-11 1.811796e-10 1.000000000
[41,] 1.793569e-10 3.587137e-10 1.000000000
[42,] 2.971402e-10 5.942803e-10 1.000000000
[43,] 5.513106e-10 1.102621e-09 0.999999999
[44,] 1.591319e-09 3.182638e-09 0.999999998
[45,] 7.415284e-09 1.483057e-08 0.999999993
[46,] 2.330512e-08 4.661023e-08 0.999999977
[47,] 6.319954e-08 1.263991e-07 0.999999937
[48,] 1.514261e-07 3.028523e-07 0.999999849
[49,] 6.035627e-07 1.207125e-06 0.999999396
[50,] 5.521185e-06 1.104237e-05 0.999994479
[51,] 1.875455e-05 3.750910e-05 0.999981245
[52,] 3.657663e-05 7.315325e-05 0.999963423
[53,] 5.515575e-05 1.103115e-04 0.999944844
[54,] 5.838004e-05 1.167601e-04 0.999941620
[55,] 4.833460e-05 9.666919e-05 0.999951665
[56,] 1.861655e-04 3.723310e-04 0.999813834
[57,] 3.432844e-04 6.865688e-04 0.999656716
[58,] 2.649400e-04 5.298800e-04 0.999735060
[59,] 1.609980e-04 3.219960e-04 0.999839002
[60,] 1.008115e-04 2.016229e-04 0.999899189
[61,] 6.737559e-05 1.347512e-04 0.999932624
[62,] 5.772990e-05 1.154598e-04 0.999942270
[63,] 1.023343e-04 2.046686e-04 0.999897666
[64,] 1.153059e-04 2.306118e-04 0.999884694
[65,] 8.168657e-05 1.633731e-04 0.999918313
[66,] 4.970019e-05 9.940038e-05 0.999950300
[67,] 4.022493e-05 8.044985e-05 0.999959775
[68,] 1.843415e-04 3.686831e-04 0.999815658
[69,] 2.196741e-04 4.393483e-04 0.999780326
[70,] 2.048663e-04 4.097325e-04 0.999795134
[71,] 1.466150e-04 2.932300e-04 0.999853385
[72,] 1.218346e-04 2.436691e-04 0.999878165
[73,] 2.302598e-04 4.605196e-04 0.999769740
[74,] 6.142675e-04 1.228535e-03 0.999385732
[75,] 1.127796e-03 2.255593e-03 0.998872204
[76,] 2.704675e-03 5.409350e-03 0.997295325
[77,] 1.628859e-02 3.257718e-02 0.983711412
[78,] 1.818624e-01 3.637249e-01 0.818137561
[79,] 8.092652e-01 3.814697e-01 0.190734836
[80,] 9.063847e-01 1.872307e-01 0.093615347
[81,] 9.189055e-01 1.621889e-01 0.081094470
[82,] 9.021107e-01 1.957785e-01 0.097889268
[83,] 9.581465e-01 8.370702e-02 0.041853509
[84,] 9.982399e-01 3.520208e-03 0.001760104
> postscript(file="/var/www/html/rcomp/tmp/10nmo1260036087.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/21fye1260036087.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/36b6y1260036087.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/4wc0k1260036087.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/5czx71260036087.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 = 95
Frequency = 1
1 2 3 4 5 6
956.803474 735.031501 694.259528 179.487555 -100.284418 694.943608
7 8 9 10 11 12
636.171635 220.399662 -254.372311 -269.144285 138.083742 183.311769
13 14 15 16 17 18
1375.539796 1246.767823 1010.995849 500.223876 505.451903 1660.679930
19 20 21 22 23 24
1336.907956 854.135983 466.364010 205.592037 311.820064 -5.951910
25 26 27 28 29 30
1199.276117 831.504144 901.732171 -89.039802 -127.811776 702.416251
31 32 33 34 35 36
610.644278 38.872305 -349.899669 -357.671642 -217.443615 -277.215588
37 38 39 40 41 42
574.012439 685.240465 128.468492 -516.303481 -558.075454 31.152572
43 44 45 46 47 48
6.380599 -499.391374 -561.163347 -849.935320 -536.707294 -845.479267
49 50 51 52 53 54
-187.251240 -97.023213 -527.795187 -972.567160 -1115.339133 -423.111106
55 56 57 58 59 60
-155.883079 -665.655053 -1120.427026 -1467.198999 -1156.970972 -1235.742945
61 62 63 64 65 66
-71.514919 -314.286892 -795.058865 -1169.830838 -1209.602812 -394.374785
67 68 69 70 71 72
-127.146758 -3661.328362 -3791.100335 -3963.872308 -3813.644281 -3353.416254
73 74 75 76 77 78
-2497.188228 -2758.960201 -2678.732174 -3412.504147 -2617.276121 -1754.048094
79 80 81 82 83 84
-1247.820067 -1195.592040 -1469.364013 -1193.135987 -640.907960 349.320067
85 86 87 88 89 90
2365.548094 3061.776121 3264.004147 2623.232174 3354.460201 4482.688228
91 92 93 94 95
4594.916254 4223.144281 3944.372308 3959.600335 3825.828362
> postscript(file="/var/www/html/rcomp/tmp/684vx1260036087.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 = 95
Frequency = 1
lag(myerror, k = 1) myerror
0 956.803474 NA
1 735.031501 956.803474
2 694.259528 735.031501
3 179.487555 694.259528
4 -100.284418 179.487555
5 694.943608 -100.284418
6 636.171635 694.943608
7 220.399662 636.171635
8 -254.372311 220.399662
9 -269.144285 -254.372311
10 138.083742 -269.144285
11 183.311769 138.083742
12 1375.539796 183.311769
13 1246.767823 1375.539796
14 1010.995849 1246.767823
15 500.223876 1010.995849
16 505.451903 500.223876
17 1660.679930 505.451903
18 1336.907956 1660.679930
19 854.135983 1336.907956
20 466.364010 854.135983
21 205.592037 466.364010
22 311.820064 205.592037
23 -5.951910 311.820064
24 1199.276117 -5.951910
25 831.504144 1199.276117
26 901.732171 831.504144
27 -89.039802 901.732171
28 -127.811776 -89.039802
29 702.416251 -127.811776
30 610.644278 702.416251
31 38.872305 610.644278
32 -349.899669 38.872305
33 -357.671642 -349.899669
34 -217.443615 -357.671642
35 -277.215588 -217.443615
36 574.012439 -277.215588
37 685.240465 574.012439
38 128.468492 685.240465
39 -516.303481 128.468492
40 -558.075454 -516.303481
41 31.152572 -558.075454
42 6.380599 31.152572
43 -499.391374 6.380599
44 -561.163347 -499.391374
45 -849.935320 -561.163347
46 -536.707294 -849.935320
47 -845.479267 -536.707294
48 -187.251240 -845.479267
49 -97.023213 -187.251240
50 -527.795187 -97.023213
51 -972.567160 -527.795187
52 -1115.339133 -972.567160
53 -423.111106 -1115.339133
54 -155.883079 -423.111106
55 -665.655053 -155.883079
56 -1120.427026 -665.655053
57 -1467.198999 -1120.427026
58 -1156.970972 -1467.198999
59 -1235.742945 -1156.970972
60 -71.514919 -1235.742945
61 -314.286892 -71.514919
62 -795.058865 -314.286892
63 -1169.830838 -795.058865
64 -1209.602812 -1169.830838
65 -394.374785 -1209.602812
66 -127.146758 -394.374785
67 -3661.328362 -127.146758
68 -3791.100335 -3661.328362
69 -3963.872308 -3791.100335
70 -3813.644281 -3963.872308
71 -3353.416254 -3813.644281
72 -2497.188228 -3353.416254
73 -2758.960201 -2497.188228
74 -2678.732174 -2758.960201
75 -3412.504147 -2678.732174
76 -2617.276121 -3412.504147
77 -1754.048094 -2617.276121
78 -1247.820067 -1754.048094
79 -1195.592040 -1247.820067
80 -1469.364013 -1195.592040
81 -1193.135987 -1469.364013
82 -640.907960 -1193.135987
83 349.320067 -640.907960
84 2365.548094 349.320067
85 3061.776121 2365.548094
86 3264.004147 3061.776121
87 2623.232174 3264.004147
88 3354.460201 2623.232174
89 4482.688228 3354.460201
90 4594.916254 4482.688228
91 4223.144281 4594.916254
92 3944.372308 4223.144281
93 3959.600335 3944.372308
94 3825.828362 3959.600335
95 NA 3825.828362
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 735.031501 956.803474
[2,] 694.259528 735.031501
[3,] 179.487555 694.259528
[4,] -100.284418 179.487555
[5,] 694.943608 -100.284418
[6,] 636.171635 694.943608
[7,] 220.399662 636.171635
[8,] -254.372311 220.399662
[9,] -269.144285 -254.372311
[10,] 138.083742 -269.144285
[11,] 183.311769 138.083742
[12,] 1375.539796 183.311769
[13,] 1246.767823 1375.539796
[14,] 1010.995849 1246.767823
[15,] 500.223876 1010.995849
[16,] 505.451903 500.223876
[17,] 1660.679930 505.451903
[18,] 1336.907956 1660.679930
[19,] 854.135983 1336.907956
[20,] 466.364010 854.135983
[21,] 205.592037 466.364010
[22,] 311.820064 205.592037
[23,] -5.951910 311.820064
[24,] 1199.276117 -5.951910
[25,] 831.504144 1199.276117
[26,] 901.732171 831.504144
[27,] -89.039802 901.732171
[28,] -127.811776 -89.039802
[29,] 702.416251 -127.811776
[30,] 610.644278 702.416251
[31,] 38.872305 610.644278
[32,] -349.899669 38.872305
[33,] -357.671642 -349.899669
[34,] -217.443615 -357.671642
[35,] -277.215588 -217.443615
[36,] 574.012439 -277.215588
[37,] 685.240465 574.012439
[38,] 128.468492 685.240465
[39,] -516.303481 128.468492
[40,] -558.075454 -516.303481
[41,] 31.152572 -558.075454
[42,] 6.380599 31.152572
[43,] -499.391374 6.380599
[44,] -561.163347 -499.391374
[45,] -849.935320 -561.163347
[46,] -536.707294 -849.935320
[47,] -845.479267 -536.707294
[48,] -187.251240 -845.479267
[49,] -97.023213 -187.251240
[50,] -527.795187 -97.023213
[51,] -972.567160 -527.795187
[52,] -1115.339133 -972.567160
[53,] -423.111106 -1115.339133
[54,] -155.883079 -423.111106
[55,] -665.655053 -155.883079
[56,] -1120.427026 -665.655053
[57,] -1467.198999 -1120.427026
[58,] -1156.970972 -1467.198999
[59,] -1235.742945 -1156.970972
[60,] -71.514919 -1235.742945
[61,] -314.286892 -71.514919
[62,] -795.058865 -314.286892
[63,] -1169.830838 -795.058865
[64,] -1209.602812 -1169.830838
[65,] -394.374785 -1209.602812
[66,] -127.146758 -394.374785
[67,] -3661.328362 -127.146758
[68,] -3791.100335 -3661.328362
[69,] -3963.872308 -3791.100335
[70,] -3813.644281 -3963.872308
[71,] -3353.416254 -3813.644281
[72,] -2497.188228 -3353.416254
[73,] -2758.960201 -2497.188228
[74,] -2678.732174 -2758.960201
[75,] -3412.504147 -2678.732174
[76,] -2617.276121 -3412.504147
[77,] -1754.048094 -2617.276121
[78,] -1247.820067 -1754.048094
[79,] -1195.592040 -1247.820067
[80,] -1469.364013 -1195.592040
[81,] -1193.135987 -1469.364013
[82,] -640.907960 -1193.135987
[83,] 349.320067 -640.907960
[84,] 2365.548094 349.320067
[85,] 3061.776121 2365.548094
[86,] 3264.004147 3061.776121
[87,] 2623.232174 3264.004147
[88,] 3354.460201 2623.232174
[89,] 4482.688228 3354.460201
[90,] 4594.916254 4482.688228
[91,] 4223.144281 4594.916254
[92,] 3944.372308 4223.144281
[93,] 3959.600335 3944.372308
[94,] 3825.828362 3959.600335
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 735.031501 956.803474
2 694.259528 735.031501
3 179.487555 694.259528
4 -100.284418 179.487555
5 694.943608 -100.284418
6 636.171635 694.943608
7 220.399662 636.171635
8 -254.372311 220.399662
9 -269.144285 -254.372311
10 138.083742 -269.144285
11 183.311769 138.083742
12 1375.539796 183.311769
13 1246.767823 1375.539796
14 1010.995849 1246.767823
15 500.223876 1010.995849
16 505.451903 500.223876
17 1660.679930 505.451903
18 1336.907956 1660.679930
19 854.135983 1336.907956
20 466.364010 854.135983
21 205.592037 466.364010
22 311.820064 205.592037
23 -5.951910 311.820064
24 1199.276117 -5.951910
25 831.504144 1199.276117
26 901.732171 831.504144
27 -89.039802 901.732171
28 -127.811776 -89.039802
29 702.416251 -127.811776
30 610.644278 702.416251
31 38.872305 610.644278
32 -349.899669 38.872305
33 -357.671642 -349.899669
34 -217.443615 -357.671642
35 -277.215588 -217.443615
36 574.012439 -277.215588
37 685.240465 574.012439
38 128.468492 685.240465
39 -516.303481 128.468492
40 -558.075454 -516.303481
41 31.152572 -558.075454
42 6.380599 31.152572
43 -499.391374 6.380599
44 -561.163347 -499.391374
45 -849.935320 -561.163347
46 -536.707294 -849.935320
47 -845.479267 -536.707294
48 -187.251240 -845.479267
49 -97.023213 -187.251240
50 -527.795187 -97.023213
51 -972.567160 -527.795187
52 -1115.339133 -972.567160
53 -423.111106 -1115.339133
54 -155.883079 -423.111106
55 -665.655053 -155.883079
56 -1120.427026 -665.655053
57 -1467.198999 -1120.427026
58 -1156.970972 -1467.198999
59 -1235.742945 -1156.970972
60 -71.514919 -1235.742945
61 -314.286892 -71.514919
62 -795.058865 -314.286892
63 -1169.830838 -795.058865
64 -1209.602812 -1169.830838
65 -394.374785 -1209.602812
66 -127.146758 -394.374785
67 -3661.328362 -127.146758
68 -3791.100335 -3661.328362
69 -3963.872308 -3791.100335
70 -3813.644281 -3963.872308
71 -3353.416254 -3813.644281
72 -2497.188228 -3353.416254
73 -2758.960201 -2497.188228
74 -2678.732174 -2758.960201
75 -3412.504147 -2678.732174
76 -2617.276121 -3412.504147
77 -1754.048094 -2617.276121
78 -1247.820067 -1754.048094
79 -1195.592040 -1247.820067
80 -1469.364013 -1195.592040
81 -1193.135987 -1469.364013
82 -640.907960 -1193.135987
83 349.320067 -640.907960
84 2365.548094 349.320067
85 3061.776121 2365.548094
86 3264.004147 3061.776121
87 2623.232174 3264.004147
88 3354.460201 2623.232174
89 4482.688228 3354.460201
90 4594.916254 4482.688228
91 4223.144281 4594.916254
92 3944.372308 4223.144281
93 3959.600335 3944.372308
94 3825.828362 3959.600335
> 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/7cxag1260036087.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/8xdic1260036087.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/9fi1c1260036087.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/10hlnd1260036087.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
+ plot(kp3:nmkm3,gqarr[,2], main='Goldfeld-Quandt test',ylab='2-sided p-value',xlab='breakpoint')
+ grid()
+ dev.off()
+ }
null device
1
>
> #Note: the /var/www/html/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/html/rcomp/createtable")
>
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Estimated Regression Equation', 1, TRUE)
> a<-table.row.end(a)
> myeq <- colnames(x)[1]
> myeq <- paste(myeq, '[t] = ', sep='')
> for (i in 1:k){
+ if (mysum$coefficients[i,1] > 0) myeq <- paste(myeq, '+', '')
+ myeq <- paste(myeq, mysum$coefficients[i,1], sep=' ')
+ if (rownames(mysum$coefficients)[i] != '(Intercept)') {
+ myeq <- paste(myeq, rownames(mysum$coefficients)[i], sep='')
+ if (rownames(mysum$coefficients)[i] != 't') myeq <- paste(myeq, '[t]', sep='')
+ }
+ }
> myeq <- paste(myeq, ' + e[t]')
> a<-table.row.start(a)
> a<-table.element(a, myeq)
> a<-table.row.end(a)
> a<-table.end(a)
> table.save(a,file="/var/www/html/rcomp/tmp/11nmgo1260036087.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/12y7ny1260036087.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/13e6911260036087.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/14a81p1260036087.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/156nq31260036087.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/16lld81260036087.tab")
+ }
>
> system("convert tmp/10nmo1260036087.ps tmp/10nmo1260036087.png")
> system("convert tmp/21fye1260036087.ps tmp/21fye1260036087.png")
> system("convert tmp/36b6y1260036087.ps tmp/36b6y1260036087.png")
> system("convert tmp/4wc0k1260036087.ps tmp/4wc0k1260036087.png")
> system("convert tmp/5czx71260036087.ps tmp/5czx71260036087.png")
> system("convert tmp/684vx1260036087.ps tmp/684vx1260036087.png")
> system("convert tmp/7cxag1260036087.ps tmp/7cxag1260036087.png")
> system("convert tmp/8xdic1260036087.ps tmp/8xdic1260036087.png")
> system("convert tmp/9fi1c1260036087.ps tmp/9fi1c1260036087.png")
> system("convert tmp/10hlnd1260036087.ps tmp/10hlnd1260036087.png")
>
>
> proc.time()
user system elapsed
2.865 1.622 3.542