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(216234,213587,209465,204045,200237,203666,241476,260307,243324,244460,233575,237217,235243,230354,227184,221678,217142,219452,256446,265845,248624,241114,229245,231805,219277,219313,212610,214771,211142,211457,240048,240636,230580,208795,197922,194596,194581,185686,178106,172608,167302,168053,202300,202388,182516,173476,166444,171297,169701,164182,161914,159612,151001,158114,186530,187069,174330,169362,166827,178037,186412,189226,191563,188906,186005,195309,223532,226899,214126),dim=c(1,69),dimnames=list(c('Werkl'),1:69))
> y <- array(NA,dim=c(1,69),dimnames=list(c('Werkl'),1:69))
> 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
Werkl M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 216234 1 0 0 0 0 0 0 0 0 0 0 1
2 213587 0 1 0 0 0 0 0 0 0 0 0 2
3 209465 0 0 1 0 0 0 0 0 0 0 0 3
4 204045 0 0 0 1 0 0 0 0 0 0 0 4
5 200237 0 0 0 0 1 0 0 0 0 0 0 5
6 203666 0 0 0 0 0 1 0 0 0 0 0 6
7 241476 0 0 0 0 0 0 1 0 0 0 0 7
8 260307 0 0 0 0 0 0 0 1 0 0 0 8
9 243324 0 0 0 0 0 0 0 0 1 0 0 9
10 244460 0 0 0 0 0 0 0 0 0 1 0 10
11 233575 0 0 0 0 0 0 0 0 0 0 1 11
12 237217 0 0 0 0 0 0 0 0 0 0 0 12
13 235243 1 0 0 0 0 0 0 0 0 0 0 13
14 230354 0 1 0 0 0 0 0 0 0 0 0 14
15 227184 0 0 1 0 0 0 0 0 0 0 0 15
16 221678 0 0 0 1 0 0 0 0 0 0 0 16
17 217142 0 0 0 0 1 0 0 0 0 0 0 17
18 219452 0 0 0 0 0 1 0 0 0 0 0 18
19 256446 0 0 0 0 0 0 1 0 0 0 0 19
20 265845 0 0 0 0 0 0 0 1 0 0 0 20
21 248624 0 0 0 0 0 0 0 0 1 0 0 21
22 241114 0 0 0 0 0 0 0 0 0 1 0 22
23 229245 0 0 0 0 0 0 0 0 0 0 1 23
24 231805 0 0 0 0 0 0 0 0 0 0 0 24
25 219277 1 0 0 0 0 0 0 0 0 0 0 25
26 219313 0 1 0 0 0 0 0 0 0 0 0 26
27 212610 0 0 1 0 0 0 0 0 0 0 0 27
28 214771 0 0 0 1 0 0 0 0 0 0 0 28
29 211142 0 0 0 0 1 0 0 0 0 0 0 29
30 211457 0 0 0 0 0 1 0 0 0 0 0 30
31 240048 0 0 0 0 0 0 1 0 0 0 0 31
32 240636 0 0 0 0 0 0 0 1 0 0 0 32
33 230580 0 0 0 0 0 0 0 0 1 0 0 33
34 208795 0 0 0 0 0 0 0 0 0 1 0 34
35 197922 0 0 0 0 0 0 0 0 0 0 1 35
36 194596 0 0 0 0 0 0 0 0 0 0 0 36
37 194581 1 0 0 0 0 0 0 0 0 0 0 37
38 185686 0 1 0 0 0 0 0 0 0 0 0 38
39 178106 0 0 1 0 0 0 0 0 0 0 0 39
40 172608 0 0 0 1 0 0 0 0 0 0 0 40
41 167302 0 0 0 0 1 0 0 0 0 0 0 41
42 168053 0 0 0 0 0 1 0 0 0 0 0 42
43 202300 0 0 0 0 0 0 1 0 0 0 0 43
44 202388 0 0 0 0 0 0 0 1 0 0 0 44
45 182516 0 0 0 0 0 0 0 0 1 0 0 45
46 173476 0 0 0 0 0 0 0 0 0 1 0 46
47 166444 0 0 0 0 0 0 0 0 0 0 1 47
48 171297 0 0 0 0 0 0 0 0 0 0 0 48
49 169701 1 0 0 0 0 0 0 0 0 0 0 49
50 164182 0 1 0 0 0 0 0 0 0 0 0 50
51 161914 0 0 1 0 0 0 0 0 0 0 0 51
52 159612 0 0 0 1 0 0 0 0 0 0 0 52
53 151001 0 0 0 0 1 0 0 0 0 0 0 53
54 158114 0 0 0 0 0 1 0 0 0 0 0 54
55 186530 0 0 0 0 0 0 1 0 0 0 0 55
56 187069 0 0 0 0 0 0 0 1 0 0 0 56
57 174330 0 0 0 0 0 0 0 0 1 0 0 57
58 169362 0 0 0 0 0 0 0 0 0 1 0 58
59 166827 0 0 0 0 0 0 0 0 0 0 1 59
60 178037 0 0 0 0 0 0 0 0 0 0 0 60
61 186412 1 0 0 0 0 0 0 0 0 0 0 61
62 189226 0 1 0 0 0 0 0 0 0 0 0 62
63 191563 0 0 1 0 0 0 0 0 0 0 0 63
64 188906 0 0 0 1 0 0 0 0 0 0 0 64
65 186005 0 0 0 0 1 0 0 0 0 0 0 65
66 195309 0 0 0 0 0 1 0 0 0 0 0 66
67 223532 0 0 0 0 0 0 1 0 0 0 0 67
68 226899 0 0 0 0 0 0 0 1 0 0 0 68
69 214126 0 0 0 0 0 0 0 0 1 0 0 69
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) M1 M2 M3 M4 M5
237091.6 -3807.6 -6032.5 -8658.5 -10903.8 -14743.9
M6 M7 M8 M9 M10 M11
-9915.2 23423.3 29850.3 15868.0 2934.3 -4746.2
t
-958.4
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-27317 -16741 1032 14018 31385
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 237091.6 9352.0 25.352 < 2e-16 ***
M1 -3807.6 11385.8 -0.334 0.7393
M2 -6032.5 11380.7 -0.530 0.5982
M3 -8658.5 11376.6 -0.761 0.4498
M4 -10903.8 11373.8 -0.959 0.3418
M5 -14743.9 11372.0 -1.297 0.2001
M6 -9915.2 11371.5 -0.872 0.3870
M7 23423.3 11372.0 2.060 0.0441 *
M8 29850.3 11373.8 2.624 0.0112 *
M9 15868.0 11376.6 1.395 0.1686
M10 2934.3 11879.3 0.247 0.8058
M11 -4746.2 11877.7 -0.400 0.6910
t -958.4 114.3 -8.386 1.8e-11 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 18780 on 56 degrees of freedom
Multiple R-squared: 0.6448, Adjusted R-squared: 0.5686
F-statistic: 8.47 on 12 and 56 DF, p-value: 7.571e-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,] 4.817191e-05 9.634381e-05 0.9999518281
[2,] 1.926200e-06 3.852401e-06 0.9999980738
[3,] 2.377170e-07 4.754339e-07 0.9999997623
[4,] 4.409472e-08 8.818945e-08 0.9999999559
[5,] 9.898835e-06 1.979767e-05 0.9999901012
[6,] 1.186083e-05 2.372167e-05 0.9999881392
[7,] 8.731105e-05 1.746221e-04 0.9999126889
[8,] 1.777018e-04 3.554036e-04 0.9998222982
[9,] 2.581335e-04 5.162670e-04 0.9997418665
[10,] 6.295135e-04 1.259027e-03 0.9993704865
[11,] 5.067116e-04 1.013423e-03 0.9994932884
[12,] 4.340745e-04 8.681489e-04 0.9995659255
[13,] 2.346499e-04 4.692999e-04 0.9997653501
[14,] 1.420489e-04 2.840977e-04 0.9998579511
[15,] 9.030537e-05 1.806107e-04 0.9999096946
[16,] 1.201429e-04 2.402859e-04 0.9998798571
[17,] 8.335224e-04 1.667045e-03 0.9991664776
[18,] 2.485849e-03 4.971698e-03 0.9975141510
[19,] 2.968107e-02 5.936215e-02 0.9703189259
[20,] 1.197979e-01 2.395958e-01 0.8802021086
[21,] 2.901968e-01 5.803937e-01 0.7098031632
[22,] 3.895810e-01 7.791620e-01 0.6104190117
[23,] 4.855359e-01 9.710719e-01 0.5144640582
[24,] 5.397751e-01 9.204497e-01 0.4602248608
[25,] 5.814770e-01 8.370460e-01 0.4185230087
[26,] 6.320437e-01 7.359126e-01 0.3679563162
[27,] 6.377658e-01 7.244684e-01 0.3622341981
[28,] 7.029009e-01 5.941981e-01 0.2970990677
[29,] 8.128652e-01 3.742696e-01 0.1871347869
[30,] 8.877916e-01 2.244169e-01 0.1122084274
[31,] 9.579895e-01 8.402100e-02 0.0420104984
[32,] 9.892156e-01 2.156872e-02 0.0107843577
[33,] 9.984897e-01 3.020639e-03 0.0015103195
[34,] 9.996737e-01 6.525684e-04 0.0003262842
[35,] 9.997449e-01 5.102814e-04 0.0002551407
[36,] 9.994858e-01 1.028487e-03 0.0005142436
[37,] 9.998421e-01 3.158041e-04 0.0001579021
[38,] 9.994772e-01 1.045654e-03 0.0005228270
> postscript(file="/var/www/html/rcomp/tmp/18ycz1259924510.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/2cg2w1259924510.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/3xvyx1259924510.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/42gja1259924510.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/5tmhg1259924510.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 = 69
Frequency = 1
1 2 3 4 5 6 7
-16091.680 -15555.347 -16093.013 -18309.347 -17318.847 -17760.180 -12330.347
8 9 10 11 12 13 14
1031.987 -1010.347 14017.789 11771.589 11625.789 14417.725 12712.059
15 16 17 18 19 20 21
13126.392 10824.059 11086.559 9526.225 14140.059 18070.392 15790.059
22 23 24 25 26 27 28
22172.195 18941.995 17714.195 9952.131 13171.464 10052.797 15417.464
29 30 31 32 33 34 35
16586.964 13031.631 9242.464 4361.797 9246.464 1353.600 -880.600
36 37 38 39 40 41 42
-7994.400 -3243.464 -8955.131 -12950.797 -15245.131 -15752.631 -18871.964
43 44 45 46 47 48 49
-17005.131 -22385.797 -27317.131 -22464.995 -20858.195 -19792.995 -16623.059
50 51 52 53 54 55 56
-18958.725 -17642.392 -16740.725 -20553.225 -17310.559 -21274.725 -26204.392
57 58 59 60 61 62 63
-24002.725 -15078.589 -8974.789 -1552.589 11588.347 17585.680 23507.013
64 65 66 67 68 69
24053.680 25951.180 31384.847 27227.680 25126.013 27293.680
> postscript(file="/var/www/html/rcomp/tmp/61y5d1259924510.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 = 69
Frequency = 1
lag(myerror, k = 1) myerror
0 -16091.680 NA
1 -15555.347 -16091.680
2 -16093.013 -15555.347
3 -18309.347 -16093.013
4 -17318.847 -18309.347
5 -17760.180 -17318.847
6 -12330.347 -17760.180
7 1031.987 -12330.347
8 -1010.347 1031.987
9 14017.789 -1010.347
10 11771.589 14017.789
11 11625.789 11771.589
12 14417.725 11625.789
13 12712.059 14417.725
14 13126.392 12712.059
15 10824.059 13126.392
16 11086.559 10824.059
17 9526.225 11086.559
18 14140.059 9526.225
19 18070.392 14140.059
20 15790.059 18070.392
21 22172.195 15790.059
22 18941.995 22172.195
23 17714.195 18941.995
24 9952.131 17714.195
25 13171.464 9952.131
26 10052.797 13171.464
27 15417.464 10052.797
28 16586.964 15417.464
29 13031.631 16586.964
30 9242.464 13031.631
31 4361.797 9242.464
32 9246.464 4361.797
33 1353.600 9246.464
34 -880.600 1353.600
35 -7994.400 -880.600
36 -3243.464 -7994.400
37 -8955.131 -3243.464
38 -12950.797 -8955.131
39 -15245.131 -12950.797
40 -15752.631 -15245.131
41 -18871.964 -15752.631
42 -17005.131 -18871.964
43 -22385.797 -17005.131
44 -27317.131 -22385.797
45 -22464.995 -27317.131
46 -20858.195 -22464.995
47 -19792.995 -20858.195
48 -16623.059 -19792.995
49 -18958.725 -16623.059
50 -17642.392 -18958.725
51 -16740.725 -17642.392
52 -20553.225 -16740.725
53 -17310.559 -20553.225
54 -21274.725 -17310.559
55 -26204.392 -21274.725
56 -24002.725 -26204.392
57 -15078.589 -24002.725
58 -8974.789 -15078.589
59 -1552.589 -8974.789
60 11588.347 -1552.589
61 17585.680 11588.347
62 23507.013 17585.680
63 24053.680 23507.013
64 25951.180 24053.680
65 31384.847 25951.180
66 27227.680 31384.847
67 25126.013 27227.680
68 27293.680 25126.013
69 NA 27293.680
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -15555.347 -16091.680
[2,] -16093.013 -15555.347
[3,] -18309.347 -16093.013
[4,] -17318.847 -18309.347
[5,] -17760.180 -17318.847
[6,] -12330.347 -17760.180
[7,] 1031.987 -12330.347
[8,] -1010.347 1031.987
[9,] 14017.789 -1010.347
[10,] 11771.589 14017.789
[11,] 11625.789 11771.589
[12,] 14417.725 11625.789
[13,] 12712.059 14417.725
[14,] 13126.392 12712.059
[15,] 10824.059 13126.392
[16,] 11086.559 10824.059
[17,] 9526.225 11086.559
[18,] 14140.059 9526.225
[19,] 18070.392 14140.059
[20,] 15790.059 18070.392
[21,] 22172.195 15790.059
[22,] 18941.995 22172.195
[23,] 17714.195 18941.995
[24,] 9952.131 17714.195
[25,] 13171.464 9952.131
[26,] 10052.797 13171.464
[27,] 15417.464 10052.797
[28,] 16586.964 15417.464
[29,] 13031.631 16586.964
[30,] 9242.464 13031.631
[31,] 4361.797 9242.464
[32,] 9246.464 4361.797
[33,] 1353.600 9246.464
[34,] -880.600 1353.600
[35,] -7994.400 -880.600
[36,] -3243.464 -7994.400
[37,] -8955.131 -3243.464
[38,] -12950.797 -8955.131
[39,] -15245.131 -12950.797
[40,] -15752.631 -15245.131
[41,] -18871.964 -15752.631
[42,] -17005.131 -18871.964
[43,] -22385.797 -17005.131
[44,] -27317.131 -22385.797
[45,] -22464.995 -27317.131
[46,] -20858.195 -22464.995
[47,] -19792.995 -20858.195
[48,] -16623.059 -19792.995
[49,] -18958.725 -16623.059
[50,] -17642.392 -18958.725
[51,] -16740.725 -17642.392
[52,] -20553.225 -16740.725
[53,] -17310.559 -20553.225
[54,] -21274.725 -17310.559
[55,] -26204.392 -21274.725
[56,] -24002.725 -26204.392
[57,] -15078.589 -24002.725
[58,] -8974.789 -15078.589
[59,] -1552.589 -8974.789
[60,] 11588.347 -1552.589
[61,] 17585.680 11588.347
[62,] 23507.013 17585.680
[63,] 24053.680 23507.013
[64,] 25951.180 24053.680
[65,] 31384.847 25951.180
[66,] 27227.680 31384.847
[67,] 25126.013 27227.680
[68,] 27293.680 25126.013
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -15555.347 -16091.680
2 -16093.013 -15555.347
3 -18309.347 -16093.013
4 -17318.847 -18309.347
5 -17760.180 -17318.847
6 -12330.347 -17760.180
7 1031.987 -12330.347
8 -1010.347 1031.987
9 14017.789 -1010.347
10 11771.589 14017.789
11 11625.789 11771.589
12 14417.725 11625.789
13 12712.059 14417.725
14 13126.392 12712.059
15 10824.059 13126.392
16 11086.559 10824.059
17 9526.225 11086.559
18 14140.059 9526.225
19 18070.392 14140.059
20 15790.059 18070.392
21 22172.195 15790.059
22 18941.995 22172.195
23 17714.195 18941.995
24 9952.131 17714.195
25 13171.464 9952.131
26 10052.797 13171.464
27 15417.464 10052.797
28 16586.964 15417.464
29 13031.631 16586.964
30 9242.464 13031.631
31 4361.797 9242.464
32 9246.464 4361.797
33 1353.600 9246.464
34 -880.600 1353.600
35 -7994.400 -880.600
36 -3243.464 -7994.400
37 -8955.131 -3243.464
38 -12950.797 -8955.131
39 -15245.131 -12950.797
40 -15752.631 -15245.131
41 -18871.964 -15752.631
42 -17005.131 -18871.964
43 -22385.797 -17005.131
44 -27317.131 -22385.797
45 -22464.995 -27317.131
46 -20858.195 -22464.995
47 -19792.995 -20858.195
48 -16623.059 -19792.995
49 -18958.725 -16623.059
50 -17642.392 -18958.725
51 -16740.725 -17642.392
52 -20553.225 -16740.725
53 -17310.559 -20553.225
54 -21274.725 -17310.559
55 -26204.392 -21274.725
56 -24002.725 -26204.392
57 -15078.589 -24002.725
58 -8974.789 -15078.589
59 -1552.589 -8974.789
60 11588.347 -1552.589
61 17585.680 11588.347
62 23507.013 17585.680
63 24053.680 23507.013
64 25951.180 24053.680
65 31384.847 25951.180
66 27227.680 31384.847
67 25126.013 27227.680
68 27293.680 25126.013
> 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/7kosr1259924510.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/8sypr1259924510.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/9wdvm1259924510.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/10ucx01259924510.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/11ku661259924510.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/12aji41259924510.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/134ewn1259924510.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/14755e1259924510.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/15kjgu1259924510.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/16hkqn1259924510.tab")
+ }
>
> system("convert tmp/18ycz1259924510.ps tmp/18ycz1259924510.png")
> system("convert tmp/2cg2w1259924510.ps tmp/2cg2w1259924510.png")
> system("convert tmp/3xvyx1259924510.ps tmp/3xvyx1259924510.png")
> system("convert tmp/42gja1259924510.ps tmp/42gja1259924510.png")
> system("convert tmp/5tmhg1259924510.ps tmp/5tmhg1259924510.png")
> system("convert tmp/61y5d1259924510.ps tmp/61y5d1259924510.png")
> system("convert tmp/7kosr1259924510.ps tmp/7kosr1259924510.png")
> system("convert tmp/8sypr1259924510.ps tmp/8sypr1259924510.png")
> system("convert tmp/9wdvm1259924510.ps tmp/9wdvm1259924510.png")
> system("convert tmp/10ucx01259924510.ps tmp/10ucx01259924510.png")
>
>
> proc.time()
user system elapsed
2.497 1.562 3.545