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(10144,112,10751,304,11752,794,13808,901,16203,1232,17432,1240,18014,1032,16956,1145,17982,1588,19435,2264,19990,2209,20154,2917,10327,243,9807,558,10862,1238,13743,1502,16458,2000,18466,2146,18810,2066,17361,2046,17411,1952,18517,2771,18525,3278,17859,4000,9499,410,9490,1107,9255,1622,10758,1986,12375,2036,14617,2400,15427,2736,14136,2901,14308,2883,15293,3747,15679,4075,16319,4996,11196,575,11169,999,12158,1411,14251,1493,16237,1846,19706,2899,18960,2372,18537,2856,19103,3468,19691,4193,19464,4440,17264,4186,8957,655,9703,1453,9166,1989,9519,2209,10535,2667,11526,3005,9630,2195,7061,2236,6021,2489,4728,2651,2657,2636,1264,2819),dim=c(2,60),dimnames=list(c('Y','X'),1:60))
> y <- array(NA,dim=c(2,60),dimnames=list(c('Y','X'),1:60))
> 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
Y X M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 10144 112 1 0 0 0 0 0 0 0 0 0 0 1
2 10751 304 0 1 0 0 0 0 0 0 0 0 0 2
3 11752 794 0 0 1 0 0 0 0 0 0 0 0 3
4 13808 901 0 0 0 1 0 0 0 0 0 0 0 4
5 16203 1232 0 0 0 0 1 0 0 0 0 0 0 5
6 17432 1240 0 0 0 0 0 1 0 0 0 0 0 6
7 18014 1032 0 0 0 0 0 0 1 0 0 0 0 7
8 16956 1145 0 0 0 0 0 0 0 1 0 0 0 8
9 17982 1588 0 0 0 0 0 0 0 0 1 0 0 9
10 19435 2264 0 0 0 0 0 0 0 0 0 1 0 10
11 19990 2209 0 0 0 0 0 0 0 0 0 0 1 11
12 20154 2917 0 0 0 0 0 0 0 0 0 0 0 12
13 10327 243 1 0 0 0 0 0 0 0 0 0 0 13
14 9807 558 0 1 0 0 0 0 0 0 0 0 0 14
15 10862 1238 0 0 1 0 0 0 0 0 0 0 0 15
16 13743 1502 0 0 0 1 0 0 0 0 0 0 0 16
17 16458 2000 0 0 0 0 1 0 0 0 0 0 0 17
18 18466 2146 0 0 0 0 0 1 0 0 0 0 0 18
19 18810 2066 0 0 0 0 0 0 1 0 0 0 0 19
20 17361 2046 0 0 0 0 0 0 0 1 0 0 0 20
21 17411 1952 0 0 0 0 0 0 0 0 1 0 0 21
22 18517 2771 0 0 0 0 0 0 0 0 0 1 0 22
23 18525 3278 0 0 0 0 0 0 0 0 0 0 1 23
24 17859 4000 0 0 0 0 0 0 0 0 0 0 0 24
25 9499 410 1 0 0 0 0 0 0 0 0 0 0 25
26 9490 1107 0 1 0 0 0 0 0 0 0 0 0 26
27 9255 1622 0 0 1 0 0 0 0 0 0 0 0 27
28 10758 1986 0 0 0 1 0 0 0 0 0 0 0 28
29 12375 2036 0 0 0 0 1 0 0 0 0 0 0 29
30 14617 2400 0 0 0 0 0 1 0 0 0 0 0 30
31 15427 2736 0 0 0 0 0 0 1 0 0 0 0 31
32 14136 2901 0 0 0 0 0 0 0 1 0 0 0 32
33 14308 2883 0 0 0 0 0 0 0 0 1 0 0 33
34 15293 3747 0 0 0 0 0 0 0 0 0 1 0 34
35 15679 4075 0 0 0 0 0 0 0 0 0 0 1 35
36 16319 4996 0 0 0 0 0 0 0 0 0 0 0 36
37 11196 575 1 0 0 0 0 0 0 0 0 0 0 37
38 11169 999 0 1 0 0 0 0 0 0 0 0 0 38
39 12158 1411 0 0 1 0 0 0 0 0 0 0 0 39
40 14251 1493 0 0 0 1 0 0 0 0 0 0 0 40
41 16237 1846 0 0 0 0 1 0 0 0 0 0 0 41
42 19706 2899 0 0 0 0 0 1 0 0 0 0 0 42
43 18960 2372 0 0 0 0 0 0 1 0 0 0 0 43
44 18537 2856 0 0 0 0 0 0 0 1 0 0 0 44
45 19103 3468 0 0 0 0 0 0 0 0 1 0 0 45
46 19691 4193 0 0 0 0 0 0 0 0 0 1 0 46
47 19464 4440 0 0 0 0 0 0 0 0 0 0 1 47
48 17264 4186 0 0 0 0 0 0 0 0 0 0 0 48
49 8957 655 1 0 0 0 0 0 0 0 0 0 0 49
50 9703 1453 0 1 0 0 0 0 0 0 0 0 0 50
51 9166 1989 0 0 1 0 0 0 0 0 0 0 0 51
52 9519 2209 0 0 0 1 0 0 0 0 0 0 0 52
53 10535 2667 0 0 0 0 1 0 0 0 0 0 0 53
54 11526 3005 0 0 0 0 0 1 0 0 0 0 0 54
55 9630 2195 0 0 0 0 0 0 1 0 0 0 0 55
56 7061 2236 0 0 0 0 0 0 0 1 0 0 0 56
57 6021 2489 0 0 0 0 0 0 0 0 1 0 0 57
58 4728 2651 0 0 0 0 0 0 0 0 0 1 0 58
59 2657 2636 0 0 0 0 0 0 0 0 0 0 1 59
60 1264 2819 0 0 0 0 0 0 0 0 0 0 0 60
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X M1 M2 M3 M4
7378.284 3.954 6461.746 4918.477 3506.721 4679.621
M5 M6 M7 M8 M9 M10
5504.753 6198.708 7252.520 5491.074 4915.842 3132.549
M11 t
2278.217 -215.720
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-4694.3 -2270.5 -314.9 2405.4 4918.9
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 7378.2836 2945.7523 2.505 0.01586 *
X 3.9538 0.8123 4.868 1.37e-05 ***
M1 6461.7459 3219.3062 2.007 0.05063 .
M2 4918.4772 2925.3443 1.681 0.09948 .
M3 3506.7208 2629.9657 1.333 0.18898
M4 4679.6206 2527.9997 1.851 0.07058 .
M5 5504.7527 2368.6591 2.324 0.02460 *
M6 6198.7080 2210.5908 2.804 0.00737 **
M7 7252.5203 2327.5742 3.116 0.00316 **
M8 5491.0736 2266.5801 2.423 0.01940 *
M9 4915.8422 2177.6977 2.257 0.02878 *
M10 3132.5486 1991.6303 1.573 0.12261
M11 2278.2174 1959.2301 1.163 0.25090
t -215.7200 28.5219 -7.563 1.31e-09 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3047 on 46 degrees of freedom
Multiple R-squared: 0.6611, Adjusted R-squared: 0.5653
F-statistic: 6.902 on 13 and 46 DF, p-value: 3.967e-07
> 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,] 3.568092e-03 7.136183e-03 0.996431908
[2,] 7.502243e-04 1.500449e-03 0.999249776
[3,] 8.567485e-05 1.713497e-04 0.999914325
[4,] 9.267555e-06 1.853511e-05 0.999990732
[5,] 1.128158e-06 2.256315e-06 0.999998872
[6,] 3.743034e-07 7.486068e-07 0.999999626
[7,] 2.740985e-06 5.481969e-06 0.999997259
[8,] 7.348182e-06 1.469636e-05 0.999992652
[9,] 1.483525e-06 2.967050e-06 0.999998516
[10,] 3.113378e-07 6.226756e-07 0.999999689
[11,] 2.257735e-07 4.515470e-07 0.999999774
[12,] 7.620231e-07 1.524046e-06 0.999999238
[13,] 2.776524e-06 5.553049e-06 0.999997223
[14,] 2.223775e-06 4.447550e-06 0.999997776
[15,] 3.168532e-06 6.337065e-06 0.999996831
[16,] 4.594691e-06 9.189381e-06 0.999995405
[17,] 4.511102e-06 9.022204e-06 0.999995489
[18,] 9.633906e-06 1.926781e-05 0.999990366
[19,] 3.006763e-05 6.013526e-05 0.999969932
[20,] 3.836778e-02 7.673556e-02 0.961632219
[21,] 3.351520e-01 6.703039e-01 0.664848035
[22,] 8.617728e-01 2.764544e-01 0.138227214
[23,] 9.929569e-01 1.408617e-02 0.007043083
[24,] 9.915180e-01 1.696405e-02 0.008482023
[25,] 9.825603e-01 3.487933e-02 0.017439665
[26,] 9.717264e-01 5.654725e-02 0.028273627
[27,] 9.188654e-01 1.622691e-01 0.081134559
> postscript(file="/var/www/html/rcomp/tmp/1kvuu1258722619.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/2gywb1258722619.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/3clu81258722619.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/4y8te1258722619.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/5y0bp1258722619.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 = 60
Frequency = 1
1 2 3 4 5 6
-3923.13624 -2316.27910 -1625.16965 -949.40711 -472.53044 246.60379
7 8 9 10 11 12
812.90404 1285.29017 1350.70366 2129.94155 3972.45228 3831.09212
13 14 15 16 17 18
-1669.44549 -1675.90701 -1682.02148 -802.00714 -665.41676 287.09167
19 20 21 22 23 24
109.30423 716.54710 1929.15664 1795.99968 869.46911 -157.24439
25 26 27 28 29 30
-569.09191 -1574.90890 -2218.64470 -3112.01137 -2302.11405 -1977.53623
31 32 33 34 35 36
-3334.10869 -3300.32070 -2266.20073 -2698.27915 -2539.07770 -3046.59942
37 38 39 40 41 42
3064.16929 3119.74247 4107.24911 4918.85690 4899.74975 3727.15238
43 44 45 46 47 48
4226.71807 3867.24063 2804.46021 2524.96140 2391.42147 3689.62667
49 50 51 52 53 54
3097.50435 2447.35254 1418.58672 -55.43128 -1459.68850 -2283.31162
55 56 57 58 59 60
-1814.81766 -2568.75720 -3818.11978 -3752.62348 -4694.26515 -4316.87498
> postscript(file="/var/www/html/rcomp/tmp/6pdu21258722619.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 = 60
Frequency = 1
lag(myerror, k = 1) myerror
0 -3923.13624 NA
1 -2316.27910 -3923.13624
2 -1625.16965 -2316.27910
3 -949.40711 -1625.16965
4 -472.53044 -949.40711
5 246.60379 -472.53044
6 812.90404 246.60379
7 1285.29017 812.90404
8 1350.70366 1285.29017
9 2129.94155 1350.70366
10 3972.45228 2129.94155
11 3831.09212 3972.45228
12 -1669.44549 3831.09212
13 -1675.90701 -1669.44549
14 -1682.02148 -1675.90701
15 -802.00714 -1682.02148
16 -665.41676 -802.00714
17 287.09167 -665.41676
18 109.30423 287.09167
19 716.54710 109.30423
20 1929.15664 716.54710
21 1795.99968 1929.15664
22 869.46911 1795.99968
23 -157.24439 869.46911
24 -569.09191 -157.24439
25 -1574.90890 -569.09191
26 -2218.64470 -1574.90890
27 -3112.01137 -2218.64470
28 -2302.11405 -3112.01137
29 -1977.53623 -2302.11405
30 -3334.10869 -1977.53623
31 -3300.32070 -3334.10869
32 -2266.20073 -3300.32070
33 -2698.27915 -2266.20073
34 -2539.07770 -2698.27915
35 -3046.59942 -2539.07770
36 3064.16929 -3046.59942
37 3119.74247 3064.16929
38 4107.24911 3119.74247
39 4918.85690 4107.24911
40 4899.74975 4918.85690
41 3727.15238 4899.74975
42 4226.71807 3727.15238
43 3867.24063 4226.71807
44 2804.46021 3867.24063
45 2524.96140 2804.46021
46 2391.42147 2524.96140
47 3689.62667 2391.42147
48 3097.50435 3689.62667
49 2447.35254 3097.50435
50 1418.58672 2447.35254
51 -55.43128 1418.58672
52 -1459.68850 -55.43128
53 -2283.31162 -1459.68850
54 -1814.81766 -2283.31162
55 -2568.75720 -1814.81766
56 -3818.11978 -2568.75720
57 -3752.62348 -3818.11978
58 -4694.26515 -3752.62348
59 -4316.87498 -4694.26515
60 NA -4316.87498
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -2316.27910 -3923.13624
[2,] -1625.16965 -2316.27910
[3,] -949.40711 -1625.16965
[4,] -472.53044 -949.40711
[5,] 246.60379 -472.53044
[6,] 812.90404 246.60379
[7,] 1285.29017 812.90404
[8,] 1350.70366 1285.29017
[9,] 2129.94155 1350.70366
[10,] 3972.45228 2129.94155
[11,] 3831.09212 3972.45228
[12,] -1669.44549 3831.09212
[13,] -1675.90701 -1669.44549
[14,] -1682.02148 -1675.90701
[15,] -802.00714 -1682.02148
[16,] -665.41676 -802.00714
[17,] 287.09167 -665.41676
[18,] 109.30423 287.09167
[19,] 716.54710 109.30423
[20,] 1929.15664 716.54710
[21,] 1795.99968 1929.15664
[22,] 869.46911 1795.99968
[23,] -157.24439 869.46911
[24,] -569.09191 -157.24439
[25,] -1574.90890 -569.09191
[26,] -2218.64470 -1574.90890
[27,] -3112.01137 -2218.64470
[28,] -2302.11405 -3112.01137
[29,] -1977.53623 -2302.11405
[30,] -3334.10869 -1977.53623
[31,] -3300.32070 -3334.10869
[32,] -2266.20073 -3300.32070
[33,] -2698.27915 -2266.20073
[34,] -2539.07770 -2698.27915
[35,] -3046.59942 -2539.07770
[36,] 3064.16929 -3046.59942
[37,] 3119.74247 3064.16929
[38,] 4107.24911 3119.74247
[39,] 4918.85690 4107.24911
[40,] 4899.74975 4918.85690
[41,] 3727.15238 4899.74975
[42,] 4226.71807 3727.15238
[43,] 3867.24063 4226.71807
[44,] 2804.46021 3867.24063
[45,] 2524.96140 2804.46021
[46,] 2391.42147 2524.96140
[47,] 3689.62667 2391.42147
[48,] 3097.50435 3689.62667
[49,] 2447.35254 3097.50435
[50,] 1418.58672 2447.35254
[51,] -55.43128 1418.58672
[52,] -1459.68850 -55.43128
[53,] -2283.31162 -1459.68850
[54,] -1814.81766 -2283.31162
[55,] -2568.75720 -1814.81766
[56,] -3818.11978 -2568.75720
[57,] -3752.62348 -3818.11978
[58,] -4694.26515 -3752.62348
[59,] -4316.87498 -4694.26515
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -2316.27910 -3923.13624
2 -1625.16965 -2316.27910
3 -949.40711 -1625.16965
4 -472.53044 -949.40711
5 246.60379 -472.53044
6 812.90404 246.60379
7 1285.29017 812.90404
8 1350.70366 1285.29017
9 2129.94155 1350.70366
10 3972.45228 2129.94155
11 3831.09212 3972.45228
12 -1669.44549 3831.09212
13 -1675.90701 -1669.44549
14 -1682.02148 -1675.90701
15 -802.00714 -1682.02148
16 -665.41676 -802.00714
17 287.09167 -665.41676
18 109.30423 287.09167
19 716.54710 109.30423
20 1929.15664 716.54710
21 1795.99968 1929.15664
22 869.46911 1795.99968
23 -157.24439 869.46911
24 -569.09191 -157.24439
25 -1574.90890 -569.09191
26 -2218.64470 -1574.90890
27 -3112.01137 -2218.64470
28 -2302.11405 -3112.01137
29 -1977.53623 -2302.11405
30 -3334.10869 -1977.53623
31 -3300.32070 -3334.10869
32 -2266.20073 -3300.32070
33 -2698.27915 -2266.20073
34 -2539.07770 -2698.27915
35 -3046.59942 -2539.07770
36 3064.16929 -3046.59942
37 3119.74247 3064.16929
38 4107.24911 3119.74247
39 4918.85690 4107.24911
40 4899.74975 4918.85690
41 3727.15238 4899.74975
42 4226.71807 3727.15238
43 3867.24063 4226.71807
44 2804.46021 3867.24063
45 2524.96140 2804.46021
46 2391.42147 2524.96140
47 3689.62667 2391.42147
48 3097.50435 3689.62667
49 2447.35254 3097.50435
50 1418.58672 2447.35254
51 -55.43128 1418.58672
52 -1459.68850 -55.43128
53 -2283.31162 -1459.68850
54 -1814.81766 -2283.31162
55 -2568.75720 -1814.81766
56 -3818.11978 -2568.75720
57 -3752.62348 -3818.11978
58 -4694.26515 -3752.62348
59 -4316.87498 -4694.26515
> 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/7napx1258722619.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/8k3ib1258722619.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/9cu391258722619.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/108tw81258722619.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/11o1pm1258722619.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/122eif1258722619.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/13awfk1258722619.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/14oobz1258722619.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/15wsy41258722619.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/16evgq1258722619.tab")
+ }
>
> system("convert tmp/1kvuu1258722619.ps tmp/1kvuu1258722619.png")
> system("convert tmp/2gywb1258722619.ps tmp/2gywb1258722619.png")
> system("convert tmp/3clu81258722619.ps tmp/3clu81258722619.png")
> system("convert tmp/4y8te1258722619.ps tmp/4y8te1258722619.png")
> system("convert tmp/5y0bp1258722619.ps tmp/5y0bp1258722619.png")
> system("convert tmp/6pdu21258722619.ps tmp/6pdu21258722619.png")
> system("convert tmp/7napx1258722619.ps tmp/7napx1258722619.png")
> system("convert tmp/8k3ib1258722619.ps tmp/8k3ib1258722619.png")
> system("convert tmp/9cu391258722619.ps tmp/9cu391258722619.png")
> system("convert tmp/108tw81258722619.ps tmp/108tw81258722619.png")
>
>
> proc.time()
user system elapsed
2.429 1.579 5.292