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(0.9059,0,0.8883,0,0.8924,0,0.8833,0,0.8700,0,0.8758,0,0.8858,0,0.9170,0,0.9554,0,0.9922,0,0.9778,0,0.9808,0,0.9811,0,1.0014,0,1.0183,0,1.0622,0,1.0773,0,1.0807,0,1.0848,0,1.1582,0,1.1663,0,1.1372,0,1.1139,0,1.1222,0,1.1692,0,1.1702,0,1.2286,0,1.2613,0,1.2646,0,1.2262,0,1.1985,0,1.2007,0,1.2138,0,1.2266,0,1.2176,0,1.2218,0,1.2490,0,1.2991,0,1.3408,0,1.3119,0,1.3014,0,1.3201,0,1.2938,0,1.2694,0,1.2165,1,1.2037,1,1.2292,1,1.2256,1,1.2015,1,1.1786,1,1.1856,1,1.2103,1,1.1938,1,1.2020,1,1.2271,1,1.2770,1,1.2650,1,1.2684,1,1.2811,1,1.2727,1,1.2611,1,1.2881,1,1.3213,1,1.2999,1,1.3074,1,1.3242,1,1.3516,1,1.3511,1,1.3419,1,1.3716,1,1.3622,1,1.3896,1,1.4227,1,1.4684,1),dim=c(2,74),dimnames=list(c('y','x'),1:74))
> y <- array(NA,dim=c(2,74),dimnames=list(c('y','x'),1:74))
> 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 0.9059 0 1 0 0 0 0 0 0 0 0 0 0 1
2 0.8883 0 0 1 0 0 0 0 0 0 0 0 0 2
3 0.8924 0 0 0 1 0 0 0 0 0 0 0 0 3
4 0.8833 0 0 0 0 1 0 0 0 0 0 0 0 4
5 0.8700 0 0 0 0 0 1 0 0 0 0 0 0 5
6 0.8758 0 0 0 0 0 0 1 0 0 0 0 0 6
7 0.8858 0 0 0 0 0 0 0 1 0 0 0 0 7
8 0.9170 0 0 0 0 0 0 0 0 1 0 0 0 8
9 0.9554 0 0 0 0 0 0 0 0 0 1 0 0 9
10 0.9922 0 0 0 0 0 0 0 0 0 0 1 0 10
11 0.9778 0 0 0 0 0 0 0 0 0 0 0 1 11
12 0.9808 0 0 0 0 0 0 0 0 0 0 0 0 12
13 0.9811 0 1 0 0 0 0 0 0 0 0 0 0 13
14 1.0014 0 0 1 0 0 0 0 0 0 0 0 0 14
15 1.0183 0 0 0 1 0 0 0 0 0 0 0 0 15
16 1.0622 0 0 0 0 1 0 0 0 0 0 0 0 16
17 1.0773 0 0 0 0 0 1 0 0 0 0 0 0 17
18 1.0807 0 0 0 0 0 0 1 0 0 0 0 0 18
19 1.0848 0 0 0 0 0 0 0 1 0 0 0 0 19
20 1.1582 0 0 0 0 0 0 0 0 1 0 0 0 20
21 1.1663 0 0 0 0 0 0 0 0 0 1 0 0 21
22 1.1372 0 0 0 0 0 0 0 0 0 0 1 0 22
23 1.1139 0 0 0 0 0 0 0 0 0 0 0 1 23
24 1.1222 0 0 0 0 0 0 0 0 0 0 0 0 24
25 1.1692 0 1 0 0 0 0 0 0 0 0 0 0 25
26 1.1702 0 0 1 0 0 0 0 0 0 0 0 0 26
27 1.2286 0 0 0 1 0 0 0 0 0 0 0 0 27
28 1.2613 0 0 0 0 1 0 0 0 0 0 0 0 28
29 1.2646 0 0 0 0 0 1 0 0 0 0 0 0 29
30 1.2262 0 0 0 0 0 0 1 0 0 0 0 0 30
31 1.1985 0 0 0 0 0 0 0 1 0 0 0 0 31
32 1.2007 0 0 0 0 0 0 0 0 1 0 0 0 32
33 1.2138 0 0 0 0 0 0 0 0 0 1 0 0 33
34 1.2266 0 0 0 0 0 0 0 0 0 0 1 0 34
35 1.2176 0 0 0 0 0 0 0 0 0 0 0 1 35
36 1.2218 0 0 0 0 0 0 0 0 0 0 0 0 36
37 1.2490 0 1 0 0 0 0 0 0 0 0 0 0 37
38 1.2991 0 0 1 0 0 0 0 0 0 0 0 0 38
39 1.3408 0 0 0 1 0 0 0 0 0 0 0 0 39
40 1.3119 0 0 0 0 1 0 0 0 0 0 0 0 40
41 1.3014 0 0 0 0 0 1 0 0 0 0 0 0 41
42 1.3201 0 0 0 0 0 0 1 0 0 0 0 0 42
43 1.2938 0 0 0 0 0 0 0 1 0 0 0 0 43
44 1.2694 0 0 0 0 0 0 0 0 1 0 0 0 44
45 1.2165 1 0 0 0 0 0 0 0 0 1 0 0 45
46 1.2037 1 0 0 0 0 0 0 0 0 0 1 0 46
47 1.2292 1 0 0 0 0 0 0 0 0 0 0 1 47
48 1.2256 1 0 0 0 0 0 0 0 0 0 0 0 48
49 1.2015 1 1 0 0 0 0 0 0 0 0 0 0 49
50 1.1786 1 0 1 0 0 0 0 0 0 0 0 0 50
51 1.1856 1 0 0 1 0 0 0 0 0 0 0 0 51
52 1.2103 1 0 0 0 1 0 0 0 0 0 0 0 52
53 1.1938 1 0 0 0 0 1 0 0 0 0 0 0 53
54 1.2020 1 0 0 0 0 0 1 0 0 0 0 0 54
55 1.2271 1 0 0 0 0 0 0 1 0 0 0 0 55
56 1.2770 1 0 0 0 0 0 0 0 1 0 0 0 56
57 1.2650 1 0 0 0 0 0 0 0 0 1 0 0 57
58 1.2684 1 0 0 0 0 0 0 0 0 0 1 0 58
59 1.2811 1 0 0 0 0 0 0 0 0 0 0 1 59
60 1.2727 1 0 0 0 0 0 0 0 0 0 0 0 60
61 1.2611 1 1 0 0 0 0 0 0 0 0 0 0 61
62 1.2881 1 0 1 0 0 0 0 0 0 0 0 0 62
63 1.3213 1 0 0 1 0 0 0 0 0 0 0 0 63
64 1.2999 1 0 0 0 1 0 0 0 0 0 0 0 64
65 1.3074 1 0 0 0 0 1 0 0 0 0 0 0 65
66 1.3242 1 0 0 0 0 0 1 0 0 0 0 0 66
67 1.3516 1 0 0 0 0 0 0 1 0 0 0 0 67
68 1.3511 1 0 0 0 0 0 0 0 1 0 0 0 68
69 1.3419 1 0 0 0 0 0 0 0 0 1 0 0 69
70 1.3716 1 0 0 0 0 0 0 0 0 0 1 0 70
71 1.3622 1 0 0 0 0 0 0 0 0 0 0 1 71
72 1.3896 1 0 0 0 0 0 0 0 0 0 0 0 72
73 1.4227 1 1 0 0 0 0 0 0 0 0 0 0 73
74 1.4684 1 0 1 0 0 0 0 0 0 0 0 0 74
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) x M1 M2 M3 M4
0.871392 -0.219909 0.004709 0.009016 0.020163 0.016654
M5 M6 M7 M8 M9 M10
0.003762 -0.004314 -0.012706 -0.001232 0.022510 0.018818
M11 t
0.005342 0.010492
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-0.062424 -0.026238 -0.006788 0.019813 0.085168
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.8713920 0.0192394 45.292 <2e-16 ***
x -0.2199086 0.0178871 -12.294 <2e-16 ***
M1 0.0047088 0.0218457 0.216 0.830
M2 0.0090164 0.0218316 0.413 0.681
M3 0.0201631 0.0227208 0.887 0.378
M4 0.0166541 0.0227038 0.734 0.466
M5 0.0037617 0.0226942 0.166 0.869
M6 -0.0043140 0.0226920 -0.190 0.850
M7 -0.0127063 0.0226971 -0.560 0.578
M8 -0.0012320 0.0227097 -0.054 0.957
M9 0.0225104 0.0226711 0.993 0.325
M10 0.0188180 0.0226526 0.831 0.409
M11 0.0053424 0.0226415 0.236 0.814
t 0.0104924 0.0004094 25.627 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.03921 on 60 degrees of freedom
Multiple R-squared: 0.9431, Adjusted R-squared: 0.9307
F-statistic: 76.44 on 13 and 60 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.8348534 0.33029313 0.16514657
[2,] 0.8407346 0.31853088 0.15926544
[3,] 0.8078570 0.38428590 0.19214295
[4,] 0.8739804 0.25203912 0.12601956
[5,] 0.8418930 0.31621391 0.15810696
[6,] 0.7825996 0.43480084 0.21740042
[7,] 0.7445675 0.51086492 0.25543246
[8,] 0.6952644 0.60947128 0.30473564
[9,] 0.6166620 0.76667595 0.38333797
[10,] 0.5426463 0.91470734 0.45735367
[11,] 0.4941081 0.98821623 0.50589188
[12,] 0.5473641 0.90527184 0.45263592
[13,] 0.6480609 0.70387822 0.35193911
[14,] 0.5879742 0.82405158 0.41202579
[15,] 0.5547967 0.89040650 0.44520325
[16,] 0.6813154 0.63736923 0.31868462
[17,] 0.7823774 0.43524528 0.21762264
[18,] 0.8111248 0.37775032 0.18887516
[19,] 0.8345542 0.33089166 0.16544583
[20,] 0.8645166 0.27096690 0.13548345
[21,] 0.8764113 0.24717744 0.12358872
[22,] 0.8286529 0.34269422 0.17134711
[23,] 0.8323259 0.33534819 0.16767409
[24,] 0.8241804 0.35163910 0.17581955
[25,] 0.8287407 0.34251857 0.17125928
[26,] 0.8595518 0.28089637 0.14044819
[27,] 0.8481396 0.30372077 0.15186039
[28,] 0.8689873 0.26202541 0.13101270
[29,] 0.8892903 0.22141936 0.11070968
[30,] 0.8681239 0.26375224 0.13187612
[31,] 0.9175419 0.16491617 0.08245809
[32,] 0.9530821 0.09383582 0.04691791
[33,] 0.9512473 0.09750536 0.04875268
[34,] 0.9447731 0.11045371 0.05522686
[35,] 0.9331086 0.13378277 0.06689138
[36,] 0.9080262 0.18394754 0.09197377
[37,] 0.8585196 0.28296073 0.14148036
[38,] 0.7799690 0.44006196 0.22003098
[39,] 0.6640885 0.67182301 0.33591150
[40,] 0.6040548 0.79189041 0.39594520
[41,] 0.5624221 0.87515583 0.43757792
> postscript(file="/var/www/html/freestat/rcomp/tmp/1qd751227376699.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/2owzn1227376699.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/3u52y1227376699.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/41km91227376699.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/5vqt81227376699.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 = 74
Frequency = 1
1 2 3 4 5
0.0193068919 -0.0130931081 -0.0306321471 -0.0467154805 -0.0576154805
6 7 8 9 10
-0.0542321471 -0.0463321471 -0.0370988138 -0.0329335811 -0.0029335811
11 12 13 14 15
-0.0143502477 -0.0165002477 -0.0314013964 -0.0259013964 -0.0306404354
16 17 18 19 20
0.0062762312 0.0237762312 0.0247595646 0.0267595646 0.0781928979
21 22 23 24 25
0.0520581306 0.0161581306 -0.0041585360 -0.0010085360 0.0307903153
26 27 28 29 30
0.0169903153 0.0537512763 0.0794679429 0.0851679429 0.0443512763
31 32 33 34 35
0.0145512763 -0.0052153904 -0.0263501577 -0.0203501577 -0.0263668243
36 37 38 39 40
-0.0273168243 -0.0153179730 0.0199820270 0.0400429880 0.0041596547
41 42 43 44 45
-0.0039403453 0.0123429880 -0.0160570120 -0.0624236787 0.0703501577
46 47 48 49 50
0.0507501577 0.0792334910 0.0704834910 0.0311823423 -0.0065176577
51 52 53 54 55
-0.0211566967 -0.0034400300 -0.0175400300 -0.0117566967 0.0112433033
56 57 58 59 60
0.0391766366 -0.0070581306 -0.0104581306 0.0052252027 -0.0083247973
61 62 63 64 65
-0.0351259459 -0.0229259459 -0.0113649850 -0.0397483183 -0.0298483183
66 67 68 69 70
-0.0154649850 0.0098350150 -0.0126316517 -0.0560664189 -0.0331664189
71 72 73 74
-0.0395830856 -0.0173330856 0.0005657658 0.0314657658
> postscript(file="/var/www/html/freestat/rcomp/tmp/6thhx1227376699.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 = 74
Frequency = 1
lag(myerror, k = 1) myerror
0 0.0193068919 NA
1 -0.0130931081 0.0193068919
2 -0.0306321471 -0.0130931081
3 -0.0467154805 -0.0306321471
4 -0.0576154805 -0.0467154805
5 -0.0542321471 -0.0576154805
6 -0.0463321471 -0.0542321471
7 -0.0370988138 -0.0463321471
8 -0.0329335811 -0.0370988138
9 -0.0029335811 -0.0329335811
10 -0.0143502477 -0.0029335811
11 -0.0165002477 -0.0143502477
12 -0.0314013964 -0.0165002477
13 -0.0259013964 -0.0314013964
14 -0.0306404354 -0.0259013964
15 0.0062762312 -0.0306404354
16 0.0237762312 0.0062762312
17 0.0247595646 0.0237762312
18 0.0267595646 0.0247595646
19 0.0781928979 0.0267595646
20 0.0520581306 0.0781928979
21 0.0161581306 0.0520581306
22 -0.0041585360 0.0161581306
23 -0.0010085360 -0.0041585360
24 0.0307903153 -0.0010085360
25 0.0169903153 0.0307903153
26 0.0537512763 0.0169903153
27 0.0794679429 0.0537512763
28 0.0851679429 0.0794679429
29 0.0443512763 0.0851679429
30 0.0145512763 0.0443512763
31 -0.0052153904 0.0145512763
32 -0.0263501577 -0.0052153904
33 -0.0203501577 -0.0263501577
34 -0.0263668243 -0.0203501577
35 -0.0273168243 -0.0263668243
36 -0.0153179730 -0.0273168243
37 0.0199820270 -0.0153179730
38 0.0400429880 0.0199820270
39 0.0041596547 0.0400429880
40 -0.0039403453 0.0041596547
41 0.0123429880 -0.0039403453
42 -0.0160570120 0.0123429880
43 -0.0624236787 -0.0160570120
44 0.0703501577 -0.0624236787
45 0.0507501577 0.0703501577
46 0.0792334910 0.0507501577
47 0.0704834910 0.0792334910
48 0.0311823423 0.0704834910
49 -0.0065176577 0.0311823423
50 -0.0211566967 -0.0065176577
51 -0.0034400300 -0.0211566967
52 -0.0175400300 -0.0034400300
53 -0.0117566967 -0.0175400300
54 0.0112433033 -0.0117566967
55 0.0391766366 0.0112433033
56 -0.0070581306 0.0391766366
57 -0.0104581306 -0.0070581306
58 0.0052252027 -0.0104581306
59 -0.0083247973 0.0052252027
60 -0.0351259459 -0.0083247973
61 -0.0229259459 -0.0351259459
62 -0.0113649850 -0.0229259459
63 -0.0397483183 -0.0113649850
64 -0.0298483183 -0.0397483183
65 -0.0154649850 -0.0298483183
66 0.0098350150 -0.0154649850
67 -0.0126316517 0.0098350150
68 -0.0560664189 -0.0126316517
69 -0.0331664189 -0.0560664189
70 -0.0395830856 -0.0331664189
71 -0.0173330856 -0.0395830856
72 0.0005657658 -0.0173330856
73 0.0314657658 0.0005657658
74 NA 0.0314657658
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -0.0130931081 0.0193068919
[2,] -0.0306321471 -0.0130931081
[3,] -0.0467154805 -0.0306321471
[4,] -0.0576154805 -0.0467154805
[5,] -0.0542321471 -0.0576154805
[6,] -0.0463321471 -0.0542321471
[7,] -0.0370988138 -0.0463321471
[8,] -0.0329335811 -0.0370988138
[9,] -0.0029335811 -0.0329335811
[10,] -0.0143502477 -0.0029335811
[11,] -0.0165002477 -0.0143502477
[12,] -0.0314013964 -0.0165002477
[13,] -0.0259013964 -0.0314013964
[14,] -0.0306404354 -0.0259013964
[15,] 0.0062762312 -0.0306404354
[16,] 0.0237762312 0.0062762312
[17,] 0.0247595646 0.0237762312
[18,] 0.0267595646 0.0247595646
[19,] 0.0781928979 0.0267595646
[20,] 0.0520581306 0.0781928979
[21,] 0.0161581306 0.0520581306
[22,] -0.0041585360 0.0161581306
[23,] -0.0010085360 -0.0041585360
[24,] 0.0307903153 -0.0010085360
[25,] 0.0169903153 0.0307903153
[26,] 0.0537512763 0.0169903153
[27,] 0.0794679429 0.0537512763
[28,] 0.0851679429 0.0794679429
[29,] 0.0443512763 0.0851679429
[30,] 0.0145512763 0.0443512763
[31,] -0.0052153904 0.0145512763
[32,] -0.0263501577 -0.0052153904
[33,] -0.0203501577 -0.0263501577
[34,] -0.0263668243 -0.0203501577
[35,] -0.0273168243 -0.0263668243
[36,] -0.0153179730 -0.0273168243
[37,] 0.0199820270 -0.0153179730
[38,] 0.0400429880 0.0199820270
[39,] 0.0041596547 0.0400429880
[40,] -0.0039403453 0.0041596547
[41,] 0.0123429880 -0.0039403453
[42,] -0.0160570120 0.0123429880
[43,] -0.0624236787 -0.0160570120
[44,] 0.0703501577 -0.0624236787
[45,] 0.0507501577 0.0703501577
[46,] 0.0792334910 0.0507501577
[47,] 0.0704834910 0.0792334910
[48,] 0.0311823423 0.0704834910
[49,] -0.0065176577 0.0311823423
[50,] -0.0211566967 -0.0065176577
[51,] -0.0034400300 -0.0211566967
[52,] -0.0175400300 -0.0034400300
[53,] -0.0117566967 -0.0175400300
[54,] 0.0112433033 -0.0117566967
[55,] 0.0391766366 0.0112433033
[56,] -0.0070581306 0.0391766366
[57,] -0.0104581306 -0.0070581306
[58,] 0.0052252027 -0.0104581306
[59,] -0.0083247973 0.0052252027
[60,] -0.0351259459 -0.0083247973
[61,] -0.0229259459 -0.0351259459
[62,] -0.0113649850 -0.0229259459
[63,] -0.0397483183 -0.0113649850
[64,] -0.0298483183 -0.0397483183
[65,] -0.0154649850 -0.0298483183
[66,] 0.0098350150 -0.0154649850
[67,] -0.0126316517 0.0098350150
[68,] -0.0560664189 -0.0126316517
[69,] -0.0331664189 -0.0560664189
[70,] -0.0395830856 -0.0331664189
[71,] -0.0173330856 -0.0395830856
[72,] 0.0005657658 -0.0173330856
[73,] 0.0314657658 0.0005657658
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -0.0130931081 0.0193068919
2 -0.0306321471 -0.0130931081
3 -0.0467154805 -0.0306321471
4 -0.0576154805 -0.0467154805
5 -0.0542321471 -0.0576154805
6 -0.0463321471 -0.0542321471
7 -0.0370988138 -0.0463321471
8 -0.0329335811 -0.0370988138
9 -0.0029335811 -0.0329335811
10 -0.0143502477 -0.0029335811
11 -0.0165002477 -0.0143502477
12 -0.0314013964 -0.0165002477
13 -0.0259013964 -0.0314013964
14 -0.0306404354 -0.0259013964
15 0.0062762312 -0.0306404354
16 0.0237762312 0.0062762312
17 0.0247595646 0.0237762312
18 0.0267595646 0.0247595646
19 0.0781928979 0.0267595646
20 0.0520581306 0.0781928979
21 0.0161581306 0.0520581306
22 -0.0041585360 0.0161581306
23 -0.0010085360 -0.0041585360
24 0.0307903153 -0.0010085360
25 0.0169903153 0.0307903153
26 0.0537512763 0.0169903153
27 0.0794679429 0.0537512763
28 0.0851679429 0.0794679429
29 0.0443512763 0.0851679429
30 0.0145512763 0.0443512763
31 -0.0052153904 0.0145512763
32 -0.0263501577 -0.0052153904
33 -0.0203501577 -0.0263501577
34 -0.0263668243 -0.0203501577
35 -0.0273168243 -0.0263668243
36 -0.0153179730 -0.0273168243
37 0.0199820270 -0.0153179730
38 0.0400429880 0.0199820270
39 0.0041596547 0.0400429880
40 -0.0039403453 0.0041596547
41 0.0123429880 -0.0039403453
42 -0.0160570120 0.0123429880
43 -0.0624236787 -0.0160570120
44 0.0703501577 -0.0624236787
45 0.0507501577 0.0703501577
46 0.0792334910 0.0507501577
47 0.0704834910 0.0792334910
48 0.0311823423 0.0704834910
49 -0.0065176577 0.0311823423
50 -0.0211566967 -0.0065176577
51 -0.0034400300 -0.0211566967
52 -0.0175400300 -0.0034400300
53 -0.0117566967 -0.0175400300
54 0.0112433033 -0.0117566967
55 0.0391766366 0.0112433033
56 -0.0070581306 0.0391766366
57 -0.0104581306 -0.0070581306
58 0.0052252027 -0.0104581306
59 -0.0083247973 0.0052252027
60 -0.0351259459 -0.0083247973
61 -0.0229259459 -0.0351259459
62 -0.0113649850 -0.0229259459
63 -0.0397483183 -0.0113649850
64 -0.0298483183 -0.0397483183
65 -0.0154649850 -0.0298483183
66 0.0098350150 -0.0154649850
67 -0.0126316517 0.0098350150
68 -0.0560664189 -0.0126316517
69 -0.0331664189 -0.0560664189
70 -0.0395830856 -0.0331664189
71 -0.0173330856 -0.0395830856
72 0.0005657658 -0.0173330856
73 0.0314657658 0.0005657658
> 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/74k621227376699.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/8wc5n1227376699.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/9k14u1227376699.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/105e121227376699.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/11pbtd1227376699.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/12cn481227376699.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/13goue1227376699.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/14vub21227376699.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/15rbcd1227376699.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/16bt2e1227376699.tab")
+ }
>
> system("convert tmp/1qd751227376699.ps tmp/1qd751227376699.png")
> system("convert tmp/2owzn1227376699.ps tmp/2owzn1227376699.png")
> system("convert tmp/3u52y1227376699.ps tmp/3u52y1227376699.png")
> system("convert tmp/41km91227376699.ps tmp/41km91227376699.png")
> system("convert tmp/5vqt81227376699.ps tmp/5vqt81227376699.png")
> system("convert tmp/6thhx1227376699.ps tmp/6thhx1227376699.png")
> system("convert tmp/74k621227376699.ps tmp/74k621227376699.png")
> system("convert tmp/8wc5n1227376699.ps tmp/8wc5n1227376699.png")
> system("convert tmp/9k14u1227376699.ps tmp/9k14u1227376699.png")
> system("convert tmp/105e121227376699.ps tmp/105e121227376699.png")
>
>
> proc.time()
user system elapsed
3.781 2.442 4.179