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.
Natural language support but running in an English locale
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(37,1,0,0,30,2,0,0,47,3,0,0,35,4,0,0,30,5,0,0,43,6,0,0,82,7,0,0,40,8,0,0,47,9,0,0,19,10,0,0,52,11,0,0,136,12,0,0,80,13,0,0,42,14,0,0,54,15,0,0,66,16,0,0,81,17,0,0,63,18,0,0,137,19,0,0,72,20,0,0,107,21,0,0,58,22,0,0,36,23,0,0,52,24,0,0,79,25,0,0,77,26,0,0,54,27,0,0,84,28,0,0,48,29,0,0,96,30,0,0,83,31,0,0,66,32,0,0,61,33,0,0,53,34,0,0,30,35,0,0,74,36,0,0,69,37,0,0,59,38,0,0,42,39,0,0,65,40,0,0,70,41,0,0,100,42,0,0,63,43,0,0,105,44,0,0,82,45,0,0,81,46,0,0,75,47,0,0,102,48,0,0,121,49,0,0,98,50,0,0,76,51,0,0,77,52,0,0,63,53,0,0,37,54,1,54,35,55,1,55,23,56,1,56,40,57,1,57,29,58,1,58,37,59,1,59,51,60,1,60,20,61,1,61,28,62,1,62,13,63,1,63,22,64,1,64,25,65,1,65,13,66,1,66,16,67,1,67,13,68,1,68,16,69,1,69,17,70,1,70,9,71,1,71,17,72,1,72,25,73,1,73,14,74,1,74,8,75,1,75,7,76,1,76,10,77,1,77,7,78,1,78,10,79,1,79,3,80,1,80),dim=c(4,80),dimnames=list(c('Y','t','D','tD'),1:80))
> y <- array(NA,dim=c(4,80),dimnames=list(c('Y','t','D','tD'),1:80))
> for (i in 1:dim(x)[1])
+ {
+ for (j in 1:dim(x)[2])
+ {
+ y[i,j] <- as.numeric(x[i,j])
+ }
+ }
> par3 = 'No Linear Trend'
> par2 = 'Do not include Seasonal Dummies'
> par1 = '1'
> #'GNU S' R Code compiled by R2WASP v. 1.0.44 ()
> #Author: Prof. Dr. P. Wessa
> #To cite this work: AUTHOR(S), (YEAR), YOUR SOFTWARE TITLE (vNUMBER) in Free Statistics Software (v$_version), Office for Research Development and Education, URL http://www.wessa.net/rwasp_YOURPAGE.wasp/
> #Source of accompanying publication: Office for Research, Development, and Education
> #Technical description: Write here your technical program description (don't use hard returns!)
> library(lattice)
> library(lmtest)
Loading required package: zoo
Attaching package: 'zoo'
The following object(s) are masked from package:base :
as.Date.numeric
> n25 <- 25 #minimum number of obs. for Goldfeld-Quandt test
> par1 <- as.numeric(par1)
> x <- t(y)
> k <- length(x[1,])
> n <- length(x[,1])
> x1 <- cbind(x[,par1], x[,1:k!=par1])
> mycolnames <- c(colnames(x)[par1], colnames(x)[1:k!=par1])
> colnames(x1) <- mycolnames #colnames(x)[par1]
> x <- x1
> if (par3 == 'First Differences'){
+ x2 <- array(0, dim=c(n-1,k), dimnames=list(1:(n-1), paste('(1-B)',colnames(x),sep='')))
+ for (i in 1:n-1) {
+ for (j in 1:k) {
+ x2[i,j] <- x[i+1,j] - x[i,j]
+ }
+ }
+ x <- x2
+ }
> if (par2 == 'Include Monthly Dummies'){
+ x2 <- array(0, dim=c(n,11), dimnames=list(1:n, paste('M', seq(1:11), sep ='')))
+ for (i in 1:11){
+ x2[seq(i,n,12),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> if (par2 == 'Include Quarterly Dummies'){
+ x2 <- array(0, dim=c(n,3), dimnames=list(1:n, paste('Q', seq(1:3), sep ='')))
+ for (i in 1:3){
+ x2[seq(i,n,4),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> k <- length(x[1,])
> if (par3 == 'Linear Trend'){
+ x <- cbind(x, c(1:n))
+ colnames(x)[k+1] <- 't'
+ }
> x
Y t D tD
1 37 1 0 0
2 30 2 0 0
3 47 3 0 0
4 35 4 0 0
5 30 5 0 0
6 43 6 0 0
7 82 7 0 0
8 40 8 0 0
9 47 9 0 0
10 19 10 0 0
11 52 11 0 0
12 136 12 0 0
13 80 13 0 0
14 42 14 0 0
15 54 15 0 0
16 66 16 0 0
17 81 17 0 0
18 63 18 0 0
19 137 19 0 0
20 72 20 0 0
21 107 21 0 0
22 58 22 0 0
23 36 23 0 0
24 52 24 0 0
25 79 25 0 0
26 77 26 0 0
27 54 27 0 0
28 84 28 0 0
29 48 29 0 0
30 96 30 0 0
31 83 31 0 0
32 66 32 0 0
33 61 33 0 0
34 53 34 0 0
35 30 35 0 0
36 74 36 0 0
37 69 37 0 0
38 59 38 0 0
39 42 39 0 0
40 65 40 0 0
41 70 41 0 0
42 100 42 0 0
43 63 43 0 0
44 105 44 0 0
45 82 45 0 0
46 81 46 0 0
47 75 47 0 0
48 102 48 0 0
49 121 49 0 0
50 98 50 0 0
51 76 51 0 0
52 77 52 0 0
53 63 53 0 0
54 37 54 1 54
55 35 55 1 55
56 23 56 1 56
57 40 57 1 57
58 29 58 1 58
59 37 59 1 59
60 51 60 1 60
61 20 61 1 61
62 28 62 1 62
63 13 63 1 63
64 22 64 1 64
65 25 65 1 65
66 13 66 1 66
67 16 67 1 67
68 13 68 1 68
69 16 69 1 69
70 17 70 1 70
71 9 71 1 71
72 17 72 1 72
73 25 73 1 73
74 14 74 1 74
75 8 75 1 75
76 7 76 1 76
77 10 77 1 77
78 7 78 1 78
79 10 79 1 79
80 3 80 1 80
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) t D tD
49.2678 0.6903 51.9473 -1.8997
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-43.428 -10.628 -2.185 7.285 78.449
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 49.2678 5.6683 8.692 5.2e-13 ***
t 0.6903 0.1827 3.779 0.000311 ***
D 51.9473 34.3719 1.511 0.134852
tD -1.8997 0.5348 -3.552 0.000660 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 20.34 on 76 degrees of freedom
Multiple R-squared: 0.6093, Adjusted R-squared: 0.5939
F-statistic: 39.51 on 3 and 76 DF, p-value: 1.710e-15
> 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.5964863 8.070274e-01 4.035137e-01
[2,] 0.5566464 8.867072e-01 4.433536e-01
[3,] 0.4347792 8.695584e-01 5.652208e-01
[4,] 0.6176159 7.647683e-01 3.823841e-01
[5,] 0.5205425 9.589151e-01 4.794575e-01
[6,] 0.9954226 9.154840e-03 4.577420e-03
[7,] 0.9918405 1.631903e-02 8.159513e-03
[8,] 0.9953391 9.321887e-03 4.660943e-03
[9,] 0.9938416 1.231675e-02 6.158376e-03
[10,] 0.9893671 2.126573e-02 1.063287e-02
[11,] 0.9833243 3.335135e-02 1.667567e-02
[12,] 0.9760684 4.786311e-02 2.393156e-02
[13,] 0.9995694 8.611979e-04 4.305990e-04
[14,] 0.9994046 1.190880e-03 5.954398e-04
[15,] 0.9998540 2.920461e-04 1.460231e-04
[16,] 0.9998971 2.058893e-04 1.029446e-04
[17,] 0.9999847 3.063532e-05 1.531766e-05
[18,] 0.9999853 2.937758e-05 1.468879e-05
[19,] 0.9999776 4.472447e-05 2.236223e-05
[20,] 0.9999667 6.660745e-05 3.330373e-05
[21,] 0.9999613 7.733633e-05 3.866817e-05
[22,] 0.9999587 8.266766e-05 4.133383e-05
[23,] 0.9999646 7.071236e-05 3.535618e-05
[24,] 0.9999887 2.262320e-05 1.131160e-05
[25,] 0.9999922 1.568214e-05 7.841071e-06
[26,] 0.9999889 2.224341e-05 1.112171e-05
[27,] 0.9999832 3.356819e-05 1.678410e-05
[28,] 0.9999774 4.524289e-05 2.262144e-05
[29,] 0.9999974 5.174519e-06 2.587260e-06
[30,] 0.9999949 1.021176e-05 5.105878e-06
[31,] 0.9999892 2.157421e-05 1.078711e-05
[32,] 0.9999826 3.473840e-05 1.736920e-05
[33,] 0.9999976 4.895122e-06 2.447561e-06
[34,] 0.9999975 4.905988e-06 2.452994e-06
[35,] 0.9999978 4.496290e-06 2.248145e-06
[36,] 0.9999970 5.950333e-06 2.975166e-06
[37,] 0.9999995 9.545321e-07 4.772661e-07
[38,] 0.9999994 1.250355e-06 6.251773e-07
[39,] 0.9999991 1.823216e-06 9.116078e-07
[40,] 0.9999994 1.273962e-06 6.369808e-07
[41,] 1.0000000 5.478248e-09 2.739124e-09
[42,] 1.0000000 1.358599e-09 6.792994e-10
[43,] 1.0000000 1.092617e-09 5.463086e-10
[44,] 1.0000000 2.450229e-09 1.225114e-09
[45,] 1.0000000 6.115074e-09 3.057537e-09
[46,] 1.0000000 2.085405e-08 1.042703e-08
[47,] 1.0000000 5.932145e-08 2.966072e-08
[48,] 0.9999999 1.978966e-07 9.894832e-08
[49,] 0.9999997 6.455441e-07 3.227720e-07
[50,] 0.9999996 8.482483e-07 4.241241e-07
[51,] 0.9999991 1.879063e-06 9.395317e-07
[52,] 0.9999972 5.691038e-06 2.845519e-06
[53,] 0.9999938 1.246793e-05 6.233966e-06
[54,] 1.0000000 5.066272e-08 2.533136e-08
[55,] 0.9999999 1.932647e-07 9.663237e-08
[56,] 0.9999998 4.178043e-07 2.089022e-07
[57,] 0.9999997 6.488542e-07 3.244271e-07
[58,] 0.9999985 2.951970e-06 1.475985e-06
[59,] 0.9999966 6.823612e-06 3.411806e-06
[60,] 0.9999897 2.066781e-05 1.033390e-05
[61,] 0.9999560 8.794405e-05 4.397203e-05
[62,] 0.9998748 2.504463e-04 1.252231e-04
[63,] 0.9995046 9.907785e-04 4.953893e-04
[64,] 0.9980076 3.984786e-03 1.992393e-03
[65,] 0.9977265 4.546994e-03 2.273497e-03
[66,] 0.9906780 1.864399e-02 9.321994e-03
[67,] 0.9932653 1.346937e-02 6.734683e-03
> postscript(file="/var/www/html/freestat/rcomp/tmp/12wal1291126405.ps",horizontal=F,onefile=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/22wal1291126405.ps",horizontal=F,onefile=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/3d5r61291126405.ps",horizontal=F,onefile=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/4d5r61291126405.ps",horizontal=F,onefile=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/5d5r61291126405.ps",horizontal=F,onefile=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 = 80
Frequency = 1
1 2 3 4 5 6
-12.95807128 -20.64836317 -4.33865506 -17.02894694 -22.71923883 -10.40953072
7 8 9 10 11 12
27.90017739 -14.79011450 -8.48040639 -37.17069827 -4.86099016 78.44871795
13 14 15 16 17 18
21.75842606 -16.93186583 -5.62215772 5.68755040 19.99725851 1.30696662
19 20 21 22 23 24
74.61667473 8.92638284 43.23609095 -6.45420094 -29.14449282 -13.83478471
25 26 27 28 29 30
12.47492340 9.78463151 -13.90566038 15.40404773 -21.28624415 26.02346396
31 32 33 34 35 36
12.33317207 -5.35711982 -11.04741171 -19.73770360 -43.42799548 -0.11828737
37 38 39 40 41 42
-5.80857926 -16.49887115 -34.18916304 -11.87945493 -7.56974682 21.73996130
43 44 45 46 47 48
-15.95033059 25.35937752 1.66908563 -0.02120626 -6.71149815 19.59820997
49 50 51 52 53 54
37.90791808 14.21762619 -8.47266570 -8.16295759 -22.85324948 1.09259259
55 56 57 58 59 60
0.30199430 -10.48860399 7.72079772 -2.06980057 7.13960114 22.34900285
61 62 63 64 65 66
-7.44159544 1.76780627 -12.02279202 -1.81339031 2.39601140 -8.39458689
67 68 69 70 71 72
-4.18518519 -5.97578348 -1.76638177 0.44301994 -6.34757835 2.86182336
73 74 75 76 77 78
12.07122507 2.28062678 -2.50997151 -2.30056980 1.90883191 0.11823362
79 80
4.32763533 -1.46296296
> postscript(file="/var/www/html/freestat/rcomp/tmp/65w8r1291126405.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> dum <- cbind(lag(myerror,k=1),myerror)
> dum
Time Series:
Start = 0
End = 80
Frequency = 1
lag(myerror, k = 1) myerror
0 -12.95807128 NA
1 -20.64836317 -12.95807128
2 -4.33865506 -20.64836317
3 -17.02894694 -4.33865506
4 -22.71923883 -17.02894694
5 -10.40953072 -22.71923883
6 27.90017739 -10.40953072
7 -14.79011450 27.90017739
8 -8.48040639 -14.79011450
9 -37.17069827 -8.48040639
10 -4.86099016 -37.17069827
11 78.44871795 -4.86099016
12 21.75842606 78.44871795
13 -16.93186583 21.75842606
14 -5.62215772 -16.93186583
15 5.68755040 -5.62215772
16 19.99725851 5.68755040
17 1.30696662 19.99725851
18 74.61667473 1.30696662
19 8.92638284 74.61667473
20 43.23609095 8.92638284
21 -6.45420094 43.23609095
22 -29.14449282 -6.45420094
23 -13.83478471 -29.14449282
24 12.47492340 -13.83478471
25 9.78463151 12.47492340
26 -13.90566038 9.78463151
27 15.40404773 -13.90566038
28 -21.28624415 15.40404773
29 26.02346396 -21.28624415
30 12.33317207 26.02346396
31 -5.35711982 12.33317207
32 -11.04741171 -5.35711982
33 -19.73770360 -11.04741171
34 -43.42799548 -19.73770360
35 -0.11828737 -43.42799548
36 -5.80857926 -0.11828737
37 -16.49887115 -5.80857926
38 -34.18916304 -16.49887115
39 -11.87945493 -34.18916304
40 -7.56974682 -11.87945493
41 21.73996130 -7.56974682
42 -15.95033059 21.73996130
43 25.35937752 -15.95033059
44 1.66908563 25.35937752
45 -0.02120626 1.66908563
46 -6.71149815 -0.02120626
47 19.59820997 -6.71149815
48 37.90791808 19.59820997
49 14.21762619 37.90791808
50 -8.47266570 14.21762619
51 -8.16295759 -8.47266570
52 -22.85324948 -8.16295759
53 1.09259259 -22.85324948
54 0.30199430 1.09259259
55 -10.48860399 0.30199430
56 7.72079772 -10.48860399
57 -2.06980057 7.72079772
58 7.13960114 -2.06980057
59 22.34900285 7.13960114
60 -7.44159544 22.34900285
61 1.76780627 -7.44159544
62 -12.02279202 1.76780627
63 -1.81339031 -12.02279202
64 2.39601140 -1.81339031
65 -8.39458689 2.39601140
66 -4.18518519 -8.39458689
67 -5.97578348 -4.18518519
68 -1.76638177 -5.97578348
69 0.44301994 -1.76638177
70 -6.34757835 0.44301994
71 2.86182336 -6.34757835
72 12.07122507 2.86182336
73 2.28062678 12.07122507
74 -2.50997151 2.28062678
75 -2.30056980 -2.50997151
76 1.90883191 -2.30056980
77 0.11823362 1.90883191
78 4.32763533 0.11823362
79 -1.46296296 4.32763533
80 NA -1.46296296
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -20.64836317 -12.95807128
[2,] -4.33865506 -20.64836317
[3,] -17.02894694 -4.33865506
[4,] -22.71923883 -17.02894694
[5,] -10.40953072 -22.71923883
[6,] 27.90017739 -10.40953072
[7,] -14.79011450 27.90017739
[8,] -8.48040639 -14.79011450
[9,] -37.17069827 -8.48040639
[10,] -4.86099016 -37.17069827
[11,] 78.44871795 -4.86099016
[12,] 21.75842606 78.44871795
[13,] -16.93186583 21.75842606
[14,] -5.62215772 -16.93186583
[15,] 5.68755040 -5.62215772
[16,] 19.99725851 5.68755040
[17,] 1.30696662 19.99725851
[18,] 74.61667473 1.30696662
[19,] 8.92638284 74.61667473
[20,] 43.23609095 8.92638284
[21,] -6.45420094 43.23609095
[22,] -29.14449282 -6.45420094
[23,] -13.83478471 -29.14449282
[24,] 12.47492340 -13.83478471
[25,] 9.78463151 12.47492340
[26,] -13.90566038 9.78463151
[27,] 15.40404773 -13.90566038
[28,] -21.28624415 15.40404773
[29,] 26.02346396 -21.28624415
[30,] 12.33317207 26.02346396
[31,] -5.35711982 12.33317207
[32,] -11.04741171 -5.35711982
[33,] -19.73770360 -11.04741171
[34,] -43.42799548 -19.73770360
[35,] -0.11828737 -43.42799548
[36,] -5.80857926 -0.11828737
[37,] -16.49887115 -5.80857926
[38,] -34.18916304 -16.49887115
[39,] -11.87945493 -34.18916304
[40,] -7.56974682 -11.87945493
[41,] 21.73996130 -7.56974682
[42,] -15.95033059 21.73996130
[43,] 25.35937752 -15.95033059
[44,] 1.66908563 25.35937752
[45,] -0.02120626 1.66908563
[46,] -6.71149815 -0.02120626
[47,] 19.59820997 -6.71149815
[48,] 37.90791808 19.59820997
[49,] 14.21762619 37.90791808
[50,] -8.47266570 14.21762619
[51,] -8.16295759 -8.47266570
[52,] -22.85324948 -8.16295759
[53,] 1.09259259 -22.85324948
[54,] 0.30199430 1.09259259
[55,] -10.48860399 0.30199430
[56,] 7.72079772 -10.48860399
[57,] -2.06980057 7.72079772
[58,] 7.13960114 -2.06980057
[59,] 22.34900285 7.13960114
[60,] -7.44159544 22.34900285
[61,] 1.76780627 -7.44159544
[62,] -12.02279202 1.76780627
[63,] -1.81339031 -12.02279202
[64,] 2.39601140 -1.81339031
[65,] -8.39458689 2.39601140
[66,] -4.18518519 -8.39458689
[67,] -5.97578348 -4.18518519
[68,] -1.76638177 -5.97578348
[69,] 0.44301994 -1.76638177
[70,] -6.34757835 0.44301994
[71,] 2.86182336 -6.34757835
[72,] 12.07122507 2.86182336
[73,] 2.28062678 12.07122507
[74,] -2.50997151 2.28062678
[75,] -2.30056980 -2.50997151
[76,] 1.90883191 -2.30056980
[77,] 0.11823362 1.90883191
[78,] 4.32763533 0.11823362
[79,] -1.46296296 4.32763533
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -20.64836317 -12.95807128
2 -4.33865506 -20.64836317
3 -17.02894694 -4.33865506
4 -22.71923883 -17.02894694
5 -10.40953072 -22.71923883
6 27.90017739 -10.40953072
7 -14.79011450 27.90017739
8 -8.48040639 -14.79011450
9 -37.17069827 -8.48040639
10 -4.86099016 -37.17069827
11 78.44871795 -4.86099016
12 21.75842606 78.44871795
13 -16.93186583 21.75842606
14 -5.62215772 -16.93186583
15 5.68755040 -5.62215772
16 19.99725851 5.68755040
17 1.30696662 19.99725851
18 74.61667473 1.30696662
19 8.92638284 74.61667473
20 43.23609095 8.92638284
21 -6.45420094 43.23609095
22 -29.14449282 -6.45420094
23 -13.83478471 -29.14449282
24 12.47492340 -13.83478471
25 9.78463151 12.47492340
26 -13.90566038 9.78463151
27 15.40404773 -13.90566038
28 -21.28624415 15.40404773
29 26.02346396 -21.28624415
30 12.33317207 26.02346396
31 -5.35711982 12.33317207
32 -11.04741171 -5.35711982
33 -19.73770360 -11.04741171
34 -43.42799548 -19.73770360
35 -0.11828737 -43.42799548
36 -5.80857926 -0.11828737
37 -16.49887115 -5.80857926
38 -34.18916304 -16.49887115
39 -11.87945493 -34.18916304
40 -7.56974682 -11.87945493
41 21.73996130 -7.56974682
42 -15.95033059 21.73996130
43 25.35937752 -15.95033059
44 1.66908563 25.35937752
45 -0.02120626 1.66908563
46 -6.71149815 -0.02120626
47 19.59820997 -6.71149815
48 37.90791808 19.59820997
49 14.21762619 37.90791808
50 -8.47266570 14.21762619
51 -8.16295759 -8.47266570
52 -22.85324948 -8.16295759
53 1.09259259 -22.85324948
54 0.30199430 1.09259259
55 -10.48860399 0.30199430
56 7.72079772 -10.48860399
57 -2.06980057 7.72079772
58 7.13960114 -2.06980057
59 22.34900285 7.13960114
60 -7.44159544 22.34900285
61 1.76780627 -7.44159544
62 -12.02279202 1.76780627
63 -1.81339031 -12.02279202
64 2.39601140 -1.81339031
65 -8.39458689 2.39601140
66 -4.18518519 -8.39458689
67 -5.97578348 -4.18518519
68 -1.76638177 -5.97578348
69 0.44301994 -1.76638177
70 -6.34757835 0.44301994
71 2.86182336 -6.34757835
72 12.07122507 2.86182336
73 2.28062678 12.07122507
74 -2.50997151 2.28062678
75 -2.30056980 -2.50997151
76 1.90883191 -2.30056980
77 0.11823362 1.90883191
78 4.32763533 0.11823362
79 -1.46296296 4.32763533
> 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/75w8r1291126405.ps",horizontal=F,onefile=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/8gopc1291126405.ps",horizontal=F,onefile=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/9gopc1291126405.ps",horizontal=F,onefile=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/109f7f1291126405.ps",horizontal=F,onefile=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/11cx531291126405.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/125p461291126405.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/134i4d1291126406.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/14qi311291126406.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/15bjjo1291126406.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/16w10c1291126406.tab")
+ }
>
> try(system("convert tmp/12wal1291126405.ps tmp/12wal1291126405.png",intern=TRUE))
character(0)
> try(system("convert tmp/22wal1291126405.ps tmp/22wal1291126405.png",intern=TRUE))
character(0)
> try(system("convert tmp/3d5r61291126405.ps tmp/3d5r61291126405.png",intern=TRUE))
character(0)
> try(system("convert tmp/4d5r61291126405.ps tmp/4d5r61291126405.png",intern=TRUE))
character(0)
> try(system("convert tmp/5d5r61291126405.ps tmp/5d5r61291126405.png",intern=TRUE))
character(0)
> try(system("convert tmp/65w8r1291126405.ps tmp/65w8r1291126405.png",intern=TRUE))
character(0)
> try(system("convert tmp/75w8r1291126405.ps tmp/75w8r1291126405.png",intern=TRUE))
character(0)
> try(system("convert tmp/8gopc1291126405.ps tmp/8gopc1291126405.png",intern=TRUE))
character(0)
> try(system("convert tmp/9gopc1291126405.ps tmp/9gopc1291126405.png",intern=TRUE))
character(0)
> try(system("convert tmp/109f7f1291126405.ps tmp/109f7f1291126405.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
4.259 2.575 5.012