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(117.3,0,118.8,0,131.3,0,125.9,0,133.1,0,147,0,145.8,0,164.4,0,149.8,0,137.7,0,151.7,0,156.8,0,180,0,180.4,0,170.4,0,191.6,0,199.5,0,218.2,1,217.5,1,205,1,194,0,199.3,0,219.3,1,211.1,1,215.2,1,240.2,1,242.2,1,240.7,1,255.4,1,253,1,218.2,1,203.7,1,205.6,1,215.6,1,188.5,1,202.9,1,214,1,230.3,1,230,1,241,1,259.6,1,247.8,1,270.3,1,289.7,1,322.7,1,315,1,320.2,1,329.5,1,360.6,1,382.2,1,435.4,1,464,1,468.8,1,403,1,351.6,1),dim=c(2,55),dimnames=list(c('Y','x'),1:55))
> y <- array(NA,dim=c(2,55),dimnames=list(c('Y','x'),1:55))
> 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 117.3 0 1 0 0 0 0 0 0 0 0 0 0 1
2 118.8 0 0 1 0 0 0 0 0 0 0 0 0 2
3 131.3 0 0 0 1 0 0 0 0 0 0 0 0 3
4 125.9 0 0 0 0 1 0 0 0 0 0 0 0 4
5 133.1 0 0 0 0 0 1 0 0 0 0 0 0 5
6 147.0 0 0 0 0 0 0 1 0 0 0 0 0 6
7 145.8 0 0 0 0 0 0 0 1 0 0 0 0 7
8 164.4 0 0 0 0 0 0 0 0 1 0 0 0 8
9 149.8 0 0 0 0 0 0 0 0 0 1 0 0 9
10 137.7 0 0 0 0 0 0 0 0 0 0 1 0 10
11 151.7 0 0 0 0 0 0 0 0 0 0 0 1 11
12 156.8 0 0 0 0 0 0 0 0 0 0 0 0 12
13 180.0 0 1 0 0 0 0 0 0 0 0 0 0 13
14 180.4 0 0 1 0 0 0 0 0 0 0 0 0 14
15 170.4 0 0 0 1 0 0 0 0 0 0 0 0 15
16 191.6 0 0 0 0 1 0 0 0 0 0 0 0 16
17 199.5 0 0 0 0 0 1 0 0 0 0 0 0 17
18 218.2 1 0 0 0 0 0 1 0 0 0 0 0 18
19 217.5 1 0 0 0 0 0 0 1 0 0 0 0 19
20 205.0 1 0 0 0 0 0 0 0 1 0 0 0 20
21 194.0 0 0 0 0 0 0 0 0 0 1 0 0 21
22 199.3 0 0 0 0 0 0 0 0 0 0 1 0 22
23 219.3 1 0 0 0 0 0 0 0 0 0 0 1 23
24 211.1 1 0 0 0 0 0 0 0 0 0 0 0 24
25 215.2 1 1 0 0 0 0 0 0 0 0 0 0 25
26 240.2 1 0 1 0 0 0 0 0 0 0 0 0 26
27 242.2 1 0 0 1 0 0 0 0 0 0 0 0 27
28 240.7 1 0 0 0 1 0 0 0 0 0 0 0 28
29 255.4 1 0 0 0 0 1 0 0 0 0 0 0 29
30 253.0 1 0 0 0 0 0 1 0 0 0 0 0 30
31 218.2 1 0 0 0 0 0 0 1 0 0 0 0 31
32 203.7 1 0 0 0 0 0 0 0 1 0 0 0 32
33 205.6 1 0 0 0 0 0 0 0 0 1 0 0 33
34 215.6 1 0 0 0 0 0 0 0 0 0 1 0 34
35 188.5 1 0 0 0 0 0 0 0 0 0 0 1 35
36 202.9 1 0 0 0 0 0 0 0 0 0 0 0 36
37 214.0 1 1 0 0 0 0 0 0 0 0 0 0 37
38 230.3 1 0 1 0 0 0 0 0 0 0 0 0 38
39 230.0 1 0 0 1 0 0 0 0 0 0 0 0 39
40 241.0 1 0 0 0 1 0 0 0 0 0 0 0 40
41 259.6 1 0 0 0 0 1 0 0 0 0 0 0 41
42 247.8 1 0 0 0 0 0 1 0 0 0 0 0 42
43 270.3 1 0 0 0 0 0 0 1 0 0 0 0 43
44 289.7 1 0 0 0 0 0 0 0 1 0 0 0 44
45 322.7 1 0 0 0 0 0 0 0 0 1 0 0 45
46 315.0 1 0 0 0 0 0 0 0 0 0 1 0 46
47 320.2 1 0 0 0 0 0 0 0 0 0 0 1 47
48 329.5 1 0 0 0 0 0 0 0 0 0 0 0 48
49 360.6 1 1 0 0 0 0 0 0 0 0 0 0 49
50 382.2 1 0 1 0 0 0 0 0 0 0 0 0 50
51 435.4 1 0 0 1 0 0 0 0 0 0 0 0 51
52 464.0 1 0 0 0 1 0 0 0 0 0 0 0 52
53 468.8 1 0 0 0 0 1 0 0 0 0 0 0 53
54 403.0 1 0 0 0 0 0 1 0 0 0 0 0 54
55 351.6 1 0 0 0 0 0 0 1 0 0 0 0 55
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) x M1 M2 M3 M4
82.6476 -34.8087 15.2127 22.5549 28.4171 33.5793
M5 M6 M7 M8 M9 M10
38.6015 30.4654 11.7276 13.0962 1.1012 -5.6416
M11 t
0.4678 5.6178
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-66.452 -17.427 4.516 18.225 90.456
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 82.6476 22.7340 3.635 0.000767 ***
x -34.8087 20.5073 -1.697 0.097204 .
M1 15.2127 27.1924 0.559 0.578901
M2 22.5549 27.1881 0.830 0.411576
M3 28.4171 27.1971 1.045 0.302211
M4 33.5793 27.2196 1.234 0.224361
M5 38.6015 27.2554 1.416 0.164244
M6 30.4654 27.1501 1.122 0.268343
M7 11.7276 27.1381 0.432 0.667899
M8 13.0962 28.7000 0.456 0.650573
M9 1.1012 28.8480 0.038 0.969735
M10 -5.6416 28.9043 -0.195 0.846215
M11 0.4678 28.6046 0.016 0.987031
t 5.6178 0.6036 9.308 1.15e-11 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 40.44 on 41 degrees of freedom
Multiple R-squared: 0.8318, Adjusted R-squared: 0.7785
F-statistic: 15.6 on 13 and 41 DF, p-value: 6.526e-12
> 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,] 8.205573e-03 1.641115e-02 0.9917944
[2,] 1.258065e-03 2.516130e-03 0.9987419
[3,] 2.039519e-04 4.079038e-04 0.9997960
[4,] 3.191514e-04 6.383028e-04 0.9996808
[5,] 8.980736e-05 1.796147e-04 0.9999102
[6,] 1.601327e-05 3.202655e-05 0.9999840
[7,] 4.431362e-06 8.862724e-06 0.9999956
[8,] 1.266345e-06 2.532690e-06 0.9999987
[9,] 1.507116e-06 3.014231e-06 0.9999985
[10,] 5.533203e-07 1.106641e-06 0.9999994
[11,] 1.696104e-07 3.392209e-07 0.9999998
[12,] 4.154261e-08 8.308521e-08 1.0000000
[13,] 1.184014e-08 2.368028e-08 1.0000000
[14,] 4.816391e-08 9.632783e-08 1.0000000
[15,] 7.119599e-05 1.423920e-04 0.9999288
[16,] 2.605940e-03 5.211880e-03 0.9973941
[17,] 5.136471e-03 1.027294e-02 0.9948635
[18,] 7.754082e-03 1.550816e-02 0.9922459
[19,] 2.122383e-02 4.244765e-02 0.9787762
[20,] 2.391271e-02 4.782542e-02 0.9760873
[21,] 1.625177e-02 3.250355e-02 0.9837482
[22,] 7.983138e-03 1.596628e-02 0.9920169
> postscript(file="/var/www/html/freestat/rcomp/tmp/1fyvs1229459144.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/2mx921229459144.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/3ny3p1229459144.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/4kxt01229459144.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/53fzk1229459144.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 = 55
Frequency = 1
1 2 3 4 5 6
13.8219215 2.3619215 3.3819215 -12.7980785 -16.2380785 0.1801759
7 8 9 10 11 12
12.1001759 23.7138227 15.4910047 4.5160047 6.7888227 6.7388227
13 14 15 16 17 18
9.1083424 -3.4516576 -24.9316576 -14.5116576 -17.2516576 38.7753248
19 20 21 22 23 24
51.1953248 31.7089716 -7.7225744 -1.2975744 41.7839716 28.4339716
25 26 27 28 29 30
11.7034912 23.7434912 14.2634912 1.9834912 6.0434912 6.1617456
31 32 33 34 35 36
-15.5182544 -37.0046076 -28.7274256 -17.6024256 -56.4296076 -47.1796076
37 38 39 40 41 42
-56.9100880 -53.5700880 -65.3500880 -65.1300880 -57.1700880 -66.4518336
43 44 45 46 47 48
-30.8318336 -18.4181867 20.9589953 14.3839953 7.8568133 12.0068133
49 50 51 52 53 54
22.2763329 30.9163329 72.6363329 90.4563329 84.6163329 21.3345873
55
-16.9454127
> postscript(file="/var/www/html/freestat/rcomp/tmp/6bhzq1229459144.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 = 55
Frequency = 1
lag(myerror, k = 1) myerror
0 13.8219215 NA
1 2.3619215 13.8219215
2 3.3819215 2.3619215
3 -12.7980785 3.3819215
4 -16.2380785 -12.7980785
5 0.1801759 -16.2380785
6 12.1001759 0.1801759
7 23.7138227 12.1001759
8 15.4910047 23.7138227
9 4.5160047 15.4910047
10 6.7888227 4.5160047
11 6.7388227 6.7888227
12 9.1083424 6.7388227
13 -3.4516576 9.1083424
14 -24.9316576 -3.4516576
15 -14.5116576 -24.9316576
16 -17.2516576 -14.5116576
17 38.7753248 -17.2516576
18 51.1953248 38.7753248
19 31.7089716 51.1953248
20 -7.7225744 31.7089716
21 -1.2975744 -7.7225744
22 41.7839716 -1.2975744
23 28.4339716 41.7839716
24 11.7034912 28.4339716
25 23.7434912 11.7034912
26 14.2634912 23.7434912
27 1.9834912 14.2634912
28 6.0434912 1.9834912
29 6.1617456 6.0434912
30 -15.5182544 6.1617456
31 -37.0046076 -15.5182544
32 -28.7274256 -37.0046076
33 -17.6024256 -28.7274256
34 -56.4296076 -17.6024256
35 -47.1796076 -56.4296076
36 -56.9100880 -47.1796076
37 -53.5700880 -56.9100880
38 -65.3500880 -53.5700880
39 -65.1300880 -65.3500880
40 -57.1700880 -65.1300880
41 -66.4518336 -57.1700880
42 -30.8318336 -66.4518336
43 -18.4181867 -30.8318336
44 20.9589953 -18.4181867
45 14.3839953 20.9589953
46 7.8568133 14.3839953
47 12.0068133 7.8568133
48 22.2763329 12.0068133
49 30.9163329 22.2763329
50 72.6363329 30.9163329
51 90.4563329 72.6363329
52 84.6163329 90.4563329
53 21.3345873 84.6163329
54 -16.9454127 21.3345873
55 NA -16.9454127
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 2.3619215 13.8219215
[2,] 3.3819215 2.3619215
[3,] -12.7980785 3.3819215
[4,] -16.2380785 -12.7980785
[5,] 0.1801759 -16.2380785
[6,] 12.1001759 0.1801759
[7,] 23.7138227 12.1001759
[8,] 15.4910047 23.7138227
[9,] 4.5160047 15.4910047
[10,] 6.7888227 4.5160047
[11,] 6.7388227 6.7888227
[12,] 9.1083424 6.7388227
[13,] -3.4516576 9.1083424
[14,] -24.9316576 -3.4516576
[15,] -14.5116576 -24.9316576
[16,] -17.2516576 -14.5116576
[17,] 38.7753248 -17.2516576
[18,] 51.1953248 38.7753248
[19,] 31.7089716 51.1953248
[20,] -7.7225744 31.7089716
[21,] -1.2975744 -7.7225744
[22,] 41.7839716 -1.2975744
[23,] 28.4339716 41.7839716
[24,] 11.7034912 28.4339716
[25,] 23.7434912 11.7034912
[26,] 14.2634912 23.7434912
[27,] 1.9834912 14.2634912
[28,] 6.0434912 1.9834912
[29,] 6.1617456 6.0434912
[30,] -15.5182544 6.1617456
[31,] -37.0046076 -15.5182544
[32,] -28.7274256 -37.0046076
[33,] -17.6024256 -28.7274256
[34,] -56.4296076 -17.6024256
[35,] -47.1796076 -56.4296076
[36,] -56.9100880 -47.1796076
[37,] -53.5700880 -56.9100880
[38,] -65.3500880 -53.5700880
[39,] -65.1300880 -65.3500880
[40,] -57.1700880 -65.1300880
[41,] -66.4518336 -57.1700880
[42,] -30.8318336 -66.4518336
[43,] -18.4181867 -30.8318336
[44,] 20.9589953 -18.4181867
[45,] 14.3839953 20.9589953
[46,] 7.8568133 14.3839953
[47,] 12.0068133 7.8568133
[48,] 22.2763329 12.0068133
[49,] 30.9163329 22.2763329
[50,] 72.6363329 30.9163329
[51,] 90.4563329 72.6363329
[52,] 84.6163329 90.4563329
[53,] 21.3345873 84.6163329
[54,] -16.9454127 21.3345873
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 2.3619215 13.8219215
2 3.3819215 2.3619215
3 -12.7980785 3.3819215
4 -16.2380785 -12.7980785
5 0.1801759 -16.2380785
6 12.1001759 0.1801759
7 23.7138227 12.1001759
8 15.4910047 23.7138227
9 4.5160047 15.4910047
10 6.7888227 4.5160047
11 6.7388227 6.7888227
12 9.1083424 6.7388227
13 -3.4516576 9.1083424
14 -24.9316576 -3.4516576
15 -14.5116576 -24.9316576
16 -17.2516576 -14.5116576
17 38.7753248 -17.2516576
18 51.1953248 38.7753248
19 31.7089716 51.1953248
20 -7.7225744 31.7089716
21 -1.2975744 -7.7225744
22 41.7839716 -1.2975744
23 28.4339716 41.7839716
24 11.7034912 28.4339716
25 23.7434912 11.7034912
26 14.2634912 23.7434912
27 1.9834912 14.2634912
28 6.0434912 1.9834912
29 6.1617456 6.0434912
30 -15.5182544 6.1617456
31 -37.0046076 -15.5182544
32 -28.7274256 -37.0046076
33 -17.6024256 -28.7274256
34 -56.4296076 -17.6024256
35 -47.1796076 -56.4296076
36 -56.9100880 -47.1796076
37 -53.5700880 -56.9100880
38 -65.3500880 -53.5700880
39 -65.1300880 -65.3500880
40 -57.1700880 -65.1300880
41 -66.4518336 -57.1700880
42 -30.8318336 -66.4518336
43 -18.4181867 -30.8318336
44 20.9589953 -18.4181867
45 14.3839953 20.9589953
46 7.8568133 14.3839953
47 12.0068133 7.8568133
48 22.2763329 12.0068133
49 30.9163329 22.2763329
50 72.6363329 30.9163329
51 90.4563329 72.6363329
52 84.6163329 90.4563329
53 21.3345873 84.6163329
54 -16.9454127 21.3345873
> 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/7m4fy1229459144.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/8hf0f1229459144.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/979gn1229459144.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/10pupl1229459144.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/11fm881229459145.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/12fqnf1229459145.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/13fv291229459145.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/145vps1229459145.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/15zpxg1229459145.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/16byqn1229459145.tab")
+ }
>
> system("convert tmp/1fyvs1229459144.ps tmp/1fyvs1229459144.png")
> system("convert tmp/2mx921229459144.ps tmp/2mx921229459144.png")
> system("convert tmp/3ny3p1229459144.ps tmp/3ny3p1229459144.png")
> system("convert tmp/4kxt01229459144.ps tmp/4kxt01229459144.png")
> system("convert tmp/53fzk1229459144.ps tmp/53fzk1229459144.png")
> system("convert tmp/6bhzq1229459144.ps tmp/6bhzq1229459144.png")
> system("convert tmp/7m4fy1229459144.ps tmp/7m4fy1229459144.png")
> system("convert tmp/8hf0f1229459144.ps tmp/8hf0f1229459144.png")
> system("convert tmp/979gn1229459144.ps tmp/979gn1229459144.png")
> system("convert tmp/10pupl1229459144.ps tmp/10pupl1229459144.png")
>
>
> proc.time()
user system elapsed
3.598 2.503 4.037