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(98.5,0,96.7,0,113.1,0,100,0,104.7,0,108.5,0,90.5,0,88.6,0,105.4,0,119.9,0,107.2,0,84.1,0,101.4,0,105.1,0,118.7,0,113.8,0,113.8,0,118.9,0,98.5,0,91,0,120.7,0,127.9,0,112.4,0,93.1,0,107.5,0,107.3,0,114.8,0,120.8,0,112.2,0,123.3,0,100.6,0,86.7,0,123.6,0,125.3,0,111.1,0,98.4,0,102.3,0,105,0,128.2,0,124.7,0,116.1,0,131.2,0,97.7,0,88.8,0,132.8,0,113.9,0,112.6,1,104.3,1,107.5,1,106,1,117.3,1,123.1,1,114.3,1,132,1,92.3,1,93.7,1,121.3,1,113.6,1,116.3,1,98.3,1,111.9,1,109.3,1,133.2,1,118,1,131.6,1,134.1,1,96.7,1,99.8,1,128.3,1,134.9,1,130.7,1,107.3,1,121.6,1,120.6,1,140.5,1,124.8,1,129.9,1,159.4,1,111,1,110.1,1,132.7,1,135,1,118.6,1,94,1,117.9,1,114.7,1,113.6,1,130.6,1,117.1,1,123.2,1,106.1,1,87.9,1),dim=c(2,92),dimnames=list(c('Y','X'),1:92))
> y <- array(NA,dim=c(2,92),dimnames=list(c('Y','X'),1:92))
> 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 X
1 98.5 0
2 96.7 0
3 113.1 0
4 100.0 0
5 104.7 0
6 108.5 0
7 90.5 0
8 88.6 0
9 105.4 0
10 119.9 0
11 107.2 0
12 84.1 0
13 101.4 0
14 105.1 0
15 118.7 0
16 113.8 0
17 113.8 0
18 118.9 0
19 98.5 0
20 91.0 0
21 120.7 0
22 127.9 0
23 112.4 0
24 93.1 0
25 107.5 0
26 107.3 0
27 114.8 0
28 120.8 0
29 112.2 0
30 123.3 0
31 100.6 0
32 86.7 0
33 123.6 0
34 125.3 0
35 111.1 0
36 98.4 0
37 102.3 0
38 105.0 0
39 128.2 0
40 124.7 0
41 116.1 0
42 131.2 0
43 97.7 0
44 88.8 0
45 132.8 0
46 113.9 0
47 112.6 1
48 104.3 1
49 107.5 1
50 106.0 1
51 117.3 1
52 123.1 1
53 114.3 1
54 132.0 1
55 92.3 1
56 93.7 1
57 121.3 1
58 113.6 1
59 116.3 1
60 98.3 1
61 111.9 1
62 109.3 1
63 133.2 1
64 118.0 1
65 131.6 1
66 134.1 1
67 96.7 1
68 99.8 1
69 128.3 1
70 134.9 1
71 130.7 1
72 107.3 1
73 121.6 1
74 120.6 1
75 140.5 1
76 124.8 1
77 129.9 1
78 159.4 1
79 111.0 1
80 110.1 1
81 132.7 1
82 135.0 1
83 118.6 1
84 94.0 1
85 117.9 1
86 114.7 1
87 113.6 1
88 130.6 1
89 117.1 1
90 123.2 1
91 106.1 1
92 87.9 1
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X
108.800 8.541
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-29.4413 -9.8913 -0.2707 10.9940 42.0587
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 108.800 2.013 54.04 < 2e-16 ***
X 8.541 2.847 3.00 0.00349 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 13.65 on 90 degrees of freedom
Multiple R-squared: 0.0909, Adjusted R-squared: 0.0808
F-statistic: 8.999 on 1 and 90 DF, p-value: 0.003493
> 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.17404693 0.34809386 0.8259531
[2,] 0.09653176 0.19306351 0.9034682
[3,] 0.12654356 0.25308712 0.8734564
[4,] 0.14925506 0.29851012 0.8507449
[5,] 0.09497679 0.18995359 0.9050232
[6,] 0.18916742 0.37833483 0.8108326
[7,] 0.12909267 0.25818534 0.8709073
[8,] 0.22420750 0.44841500 0.7757925
[9,] 0.15879619 0.31759237 0.8412038
[10,] 0.11114426 0.22228851 0.8888557
[11,] 0.14653158 0.29306316 0.8534684
[12,] 0.13024146 0.26048292 0.8697585
[13,] 0.11211495 0.22422991 0.8878850
[14,] 0.12174180 0.24348361 0.8782582
[15,] 0.09660197 0.19320395 0.9033980
[16,] 0.10992843 0.21985687 0.8900716
[17,] 0.13069886 0.26139772 0.8693011
[18,] 0.21960252 0.43920505 0.7803975
[19,] 0.17812730 0.35625460 0.8218727
[20,] 0.18569429 0.37138857 0.8143057
[21,] 0.14313720 0.28627439 0.8568628
[22,] 0.10811149 0.21622298 0.8918885
[23,] 0.08963894 0.17927788 0.9103611
[24,] 0.09249157 0.18498314 0.9075084
[25,] 0.07043815 0.14087631 0.9295618
[26,] 0.08023818 0.16047635 0.9197618
[27,] 0.06561442 0.13122884 0.9343856
[28,] 0.11319105 0.22638209 0.8868090
[29,] 0.12498357 0.24996714 0.8750164
[30,] 0.14500764 0.29001529 0.8549924
[31,] 0.11310518 0.22621035 0.8868948
[32,] 0.10405896 0.20811792 0.8959410
[33,] 0.08683319 0.17366639 0.9131668
[34,] 0.06910843 0.13821687 0.9308916
[35,] 0.09151562 0.18303125 0.9084844
[36,] 0.09809316 0.19618632 0.9019068
[37,] 0.07925635 0.15851270 0.9207437
[38,] 0.12246222 0.24492445 0.8775378
[39,] 0.11356781 0.22713563 0.8864322
[40,] 0.18093363 0.36186726 0.8190664
[41,] 0.23452687 0.46905374 0.7654731
[42,] 0.19365887 0.38731775 0.8063411
[43,] 0.15651110 0.31302220 0.8434889
[44,] 0.13923711 0.27847423 0.8607629
[45,] 0.11555008 0.23110015 0.8844499
[46,] 0.09762230 0.19524461 0.9023777
[47,] 0.07872142 0.15744284 0.9212786
[48,] 0.06807258 0.13614516 0.9319274
[49,] 0.05100614 0.10201227 0.9489939
[50,] 0.05819115 0.11638231 0.9418088
[51,] 0.09984435 0.19968870 0.9001556
[52,] 0.14663856 0.29327712 0.8533614
[53,] 0.12239308 0.24478615 0.8776069
[54,] 0.09633702 0.19267403 0.9036630
[55,] 0.07397824 0.14795649 0.9260218
[56,] 0.09126590 0.18253180 0.9087341
[57,] 0.07240538 0.14481075 0.9275946
[58,] 0.05984280 0.11968561 0.9401572
[59,] 0.06983829 0.13967657 0.9301617
[60,] 0.05207699 0.10415399 0.9479230
[61,] 0.05346754 0.10693508 0.9465325
[62,] 0.06063287 0.12126573 0.9393671
[63,] 0.08902352 0.17804704 0.9109765
[64,] 0.11146412 0.22292824 0.8885359
[65,] 0.09688281 0.19376563 0.9031172
[66,] 0.10717037 0.21434074 0.8928296
[67,] 0.09865747 0.19731493 0.9013425
[68,] 0.08762044 0.17524088 0.9123796
[69,] 0.06333780 0.12667560 0.9366622
[70,] 0.04398122 0.08796245 0.9560188
[71,] 0.06808850 0.13617699 0.9319115
[72,] 0.04964828 0.09929655 0.9503517
[73,] 0.04210178 0.08420355 0.9578982
[74,] 0.42568743 0.85137485 0.5743126
[75,] 0.34895558 0.69791116 0.6510444
[76,] 0.28022749 0.56045497 0.7197725
[77,] 0.31279502 0.62559005 0.6872050
[78,] 0.42803305 0.85606610 0.5719670
[79,] 0.34268930 0.68537860 0.6573107
[80,] 0.43542009 0.87084019 0.5645799
[81,] 0.32269409 0.64538818 0.6773059
[82,] 0.21034214 0.42068427 0.7896579
[83,] 0.11759115 0.23518231 0.8824088
> postscript(file="/var/www/html/freestat/rcomp/tmp/11iwd1227731996.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/2l53j1227731996.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/3rkht1227731996.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/4yowk1227731996.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/581i71227731996.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 = 92
Frequency = 1
1 2 3 4 5 6
-10.30000000 -12.10000000 4.30000000 -8.80000000 -4.10000000 -0.30000000
7 8 9 10 11 12
-18.30000000 -20.20000000 -3.40000000 11.10000000 -1.60000000 -24.70000000
13 14 15 16 17 18
-7.40000000 -3.70000000 9.90000000 5.00000000 5.00000000 10.10000000
19 20 21 22 23 24
-10.30000000 -17.80000000 11.90000000 19.10000000 3.60000000 -15.70000000
25 26 27 28 29 30
-1.30000000 -1.50000000 6.00000000 12.00000000 3.40000000 14.50000000
31 32 33 34 35 36
-8.20000000 -22.10000000 14.80000000 16.50000000 2.30000000 -10.40000000
37 38 39 40 41 42
-6.50000000 -3.80000000 19.40000000 15.90000000 7.30000000 22.40000000
43 44 45 46 47 48
-11.10000000 -20.00000000 24.00000000 5.10000000 -4.74130435 -13.04130435
49 50 51 52 53 54
-9.84130435 -11.34130435 -0.04130435 5.75869565 -3.04130435 14.65869565
55 56 57 58 59 60
-25.04130435 -23.64130435 3.95869565 -3.74130435 -1.04130435 -19.04130435
61 62 63 64 65 66
-5.44130435 -8.04130435 15.85869565 0.65869565 14.25869565 16.75869565
67 68 69 70 71 72
-20.64130435 -17.54130435 10.95869565 17.55869565 13.35869565 -10.04130435
73 74 75 76 77 78
4.25869565 3.25869565 23.15869565 7.45869565 12.55869565 42.05869565
79 80 81 82 83 84
-6.34130435 -7.24130435 15.35869565 17.65869565 1.25869565 -23.34130435
85 86 87 88 89 90
0.55869565 -2.64130435 -3.74130435 13.25869565 -0.24130435 5.85869565
91 92
-11.24130435 -29.44130435
> postscript(file="/var/www/html/freestat/rcomp/tmp/6dkd71227731996.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 = 92
Frequency = 1
lag(myerror, k = 1) myerror
0 -10.30000000 NA
1 -12.10000000 -10.30000000
2 4.30000000 -12.10000000
3 -8.80000000 4.30000000
4 -4.10000000 -8.80000000
5 -0.30000000 -4.10000000
6 -18.30000000 -0.30000000
7 -20.20000000 -18.30000000
8 -3.40000000 -20.20000000
9 11.10000000 -3.40000000
10 -1.60000000 11.10000000
11 -24.70000000 -1.60000000
12 -7.40000000 -24.70000000
13 -3.70000000 -7.40000000
14 9.90000000 -3.70000000
15 5.00000000 9.90000000
16 5.00000000 5.00000000
17 10.10000000 5.00000000
18 -10.30000000 10.10000000
19 -17.80000000 -10.30000000
20 11.90000000 -17.80000000
21 19.10000000 11.90000000
22 3.60000000 19.10000000
23 -15.70000000 3.60000000
24 -1.30000000 -15.70000000
25 -1.50000000 -1.30000000
26 6.00000000 -1.50000000
27 12.00000000 6.00000000
28 3.40000000 12.00000000
29 14.50000000 3.40000000
30 -8.20000000 14.50000000
31 -22.10000000 -8.20000000
32 14.80000000 -22.10000000
33 16.50000000 14.80000000
34 2.30000000 16.50000000
35 -10.40000000 2.30000000
36 -6.50000000 -10.40000000
37 -3.80000000 -6.50000000
38 19.40000000 -3.80000000
39 15.90000000 19.40000000
40 7.30000000 15.90000000
41 22.40000000 7.30000000
42 -11.10000000 22.40000000
43 -20.00000000 -11.10000000
44 24.00000000 -20.00000000
45 5.10000000 24.00000000
46 -4.74130435 5.10000000
47 -13.04130435 -4.74130435
48 -9.84130435 -13.04130435
49 -11.34130435 -9.84130435
50 -0.04130435 -11.34130435
51 5.75869565 -0.04130435
52 -3.04130435 5.75869565
53 14.65869565 -3.04130435
54 -25.04130435 14.65869565
55 -23.64130435 -25.04130435
56 3.95869565 -23.64130435
57 -3.74130435 3.95869565
58 -1.04130435 -3.74130435
59 -19.04130435 -1.04130435
60 -5.44130435 -19.04130435
61 -8.04130435 -5.44130435
62 15.85869565 -8.04130435
63 0.65869565 15.85869565
64 14.25869565 0.65869565
65 16.75869565 14.25869565
66 -20.64130435 16.75869565
67 -17.54130435 -20.64130435
68 10.95869565 -17.54130435
69 17.55869565 10.95869565
70 13.35869565 17.55869565
71 -10.04130435 13.35869565
72 4.25869565 -10.04130435
73 3.25869565 4.25869565
74 23.15869565 3.25869565
75 7.45869565 23.15869565
76 12.55869565 7.45869565
77 42.05869565 12.55869565
78 -6.34130435 42.05869565
79 -7.24130435 -6.34130435
80 15.35869565 -7.24130435
81 17.65869565 15.35869565
82 1.25869565 17.65869565
83 -23.34130435 1.25869565
84 0.55869565 -23.34130435
85 -2.64130435 0.55869565
86 -3.74130435 -2.64130435
87 13.25869565 -3.74130435
88 -0.24130435 13.25869565
89 5.85869565 -0.24130435
90 -11.24130435 5.85869565
91 -29.44130435 -11.24130435
92 NA -29.44130435
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -12.10000000 -10.30000000
[2,] 4.30000000 -12.10000000
[3,] -8.80000000 4.30000000
[4,] -4.10000000 -8.80000000
[5,] -0.30000000 -4.10000000
[6,] -18.30000000 -0.30000000
[7,] -20.20000000 -18.30000000
[8,] -3.40000000 -20.20000000
[9,] 11.10000000 -3.40000000
[10,] -1.60000000 11.10000000
[11,] -24.70000000 -1.60000000
[12,] -7.40000000 -24.70000000
[13,] -3.70000000 -7.40000000
[14,] 9.90000000 -3.70000000
[15,] 5.00000000 9.90000000
[16,] 5.00000000 5.00000000
[17,] 10.10000000 5.00000000
[18,] -10.30000000 10.10000000
[19,] -17.80000000 -10.30000000
[20,] 11.90000000 -17.80000000
[21,] 19.10000000 11.90000000
[22,] 3.60000000 19.10000000
[23,] -15.70000000 3.60000000
[24,] -1.30000000 -15.70000000
[25,] -1.50000000 -1.30000000
[26,] 6.00000000 -1.50000000
[27,] 12.00000000 6.00000000
[28,] 3.40000000 12.00000000
[29,] 14.50000000 3.40000000
[30,] -8.20000000 14.50000000
[31,] -22.10000000 -8.20000000
[32,] 14.80000000 -22.10000000
[33,] 16.50000000 14.80000000
[34,] 2.30000000 16.50000000
[35,] -10.40000000 2.30000000
[36,] -6.50000000 -10.40000000
[37,] -3.80000000 -6.50000000
[38,] 19.40000000 -3.80000000
[39,] 15.90000000 19.40000000
[40,] 7.30000000 15.90000000
[41,] 22.40000000 7.30000000
[42,] -11.10000000 22.40000000
[43,] -20.00000000 -11.10000000
[44,] 24.00000000 -20.00000000
[45,] 5.10000000 24.00000000
[46,] -4.74130435 5.10000000
[47,] -13.04130435 -4.74130435
[48,] -9.84130435 -13.04130435
[49,] -11.34130435 -9.84130435
[50,] -0.04130435 -11.34130435
[51,] 5.75869565 -0.04130435
[52,] -3.04130435 5.75869565
[53,] 14.65869565 -3.04130435
[54,] -25.04130435 14.65869565
[55,] -23.64130435 -25.04130435
[56,] 3.95869565 -23.64130435
[57,] -3.74130435 3.95869565
[58,] -1.04130435 -3.74130435
[59,] -19.04130435 -1.04130435
[60,] -5.44130435 -19.04130435
[61,] -8.04130435 -5.44130435
[62,] 15.85869565 -8.04130435
[63,] 0.65869565 15.85869565
[64,] 14.25869565 0.65869565
[65,] 16.75869565 14.25869565
[66,] -20.64130435 16.75869565
[67,] -17.54130435 -20.64130435
[68,] 10.95869565 -17.54130435
[69,] 17.55869565 10.95869565
[70,] 13.35869565 17.55869565
[71,] -10.04130435 13.35869565
[72,] 4.25869565 -10.04130435
[73,] 3.25869565 4.25869565
[74,] 23.15869565 3.25869565
[75,] 7.45869565 23.15869565
[76,] 12.55869565 7.45869565
[77,] 42.05869565 12.55869565
[78,] -6.34130435 42.05869565
[79,] -7.24130435 -6.34130435
[80,] 15.35869565 -7.24130435
[81,] 17.65869565 15.35869565
[82,] 1.25869565 17.65869565
[83,] -23.34130435 1.25869565
[84,] 0.55869565 -23.34130435
[85,] -2.64130435 0.55869565
[86,] -3.74130435 -2.64130435
[87,] 13.25869565 -3.74130435
[88,] -0.24130435 13.25869565
[89,] 5.85869565 -0.24130435
[90,] -11.24130435 5.85869565
[91,] -29.44130435 -11.24130435
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -12.10000000 -10.30000000
2 4.30000000 -12.10000000
3 -8.80000000 4.30000000
4 -4.10000000 -8.80000000
5 -0.30000000 -4.10000000
6 -18.30000000 -0.30000000
7 -20.20000000 -18.30000000
8 -3.40000000 -20.20000000
9 11.10000000 -3.40000000
10 -1.60000000 11.10000000
11 -24.70000000 -1.60000000
12 -7.40000000 -24.70000000
13 -3.70000000 -7.40000000
14 9.90000000 -3.70000000
15 5.00000000 9.90000000
16 5.00000000 5.00000000
17 10.10000000 5.00000000
18 -10.30000000 10.10000000
19 -17.80000000 -10.30000000
20 11.90000000 -17.80000000
21 19.10000000 11.90000000
22 3.60000000 19.10000000
23 -15.70000000 3.60000000
24 -1.30000000 -15.70000000
25 -1.50000000 -1.30000000
26 6.00000000 -1.50000000
27 12.00000000 6.00000000
28 3.40000000 12.00000000
29 14.50000000 3.40000000
30 -8.20000000 14.50000000
31 -22.10000000 -8.20000000
32 14.80000000 -22.10000000
33 16.50000000 14.80000000
34 2.30000000 16.50000000
35 -10.40000000 2.30000000
36 -6.50000000 -10.40000000
37 -3.80000000 -6.50000000
38 19.40000000 -3.80000000
39 15.90000000 19.40000000
40 7.30000000 15.90000000
41 22.40000000 7.30000000
42 -11.10000000 22.40000000
43 -20.00000000 -11.10000000
44 24.00000000 -20.00000000
45 5.10000000 24.00000000
46 -4.74130435 5.10000000
47 -13.04130435 -4.74130435
48 -9.84130435 -13.04130435
49 -11.34130435 -9.84130435
50 -0.04130435 -11.34130435
51 5.75869565 -0.04130435
52 -3.04130435 5.75869565
53 14.65869565 -3.04130435
54 -25.04130435 14.65869565
55 -23.64130435 -25.04130435
56 3.95869565 -23.64130435
57 -3.74130435 3.95869565
58 -1.04130435 -3.74130435
59 -19.04130435 -1.04130435
60 -5.44130435 -19.04130435
61 -8.04130435 -5.44130435
62 15.85869565 -8.04130435
63 0.65869565 15.85869565
64 14.25869565 0.65869565
65 16.75869565 14.25869565
66 -20.64130435 16.75869565
67 -17.54130435 -20.64130435
68 10.95869565 -17.54130435
69 17.55869565 10.95869565
70 13.35869565 17.55869565
71 -10.04130435 13.35869565
72 4.25869565 -10.04130435
73 3.25869565 4.25869565
74 23.15869565 3.25869565
75 7.45869565 23.15869565
76 12.55869565 7.45869565
77 42.05869565 12.55869565
78 -6.34130435 42.05869565
79 -7.24130435 -6.34130435
80 15.35869565 -7.24130435
81 17.65869565 15.35869565
82 1.25869565 17.65869565
83 -23.34130435 1.25869565
84 0.55869565 -23.34130435
85 -2.64130435 0.55869565
86 -3.74130435 -2.64130435
87 13.25869565 -3.74130435
88 -0.24130435 13.25869565
89 5.85869565 -0.24130435
90 -11.24130435 5.85869565
91 -29.44130435 -11.24130435
> 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/76lz21227731996.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/85q0q1227731996.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/9qm981227731996.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')
hat values (leverages) are all = 0.02173913
and there are no factor predictors; no plot no. 5
> par(opar)
> dev.off()
null device
1
> if (n > n25) {
+ postscript(file="/var/www/html/freestat/rcomp/tmp/10nc731227731996.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/11ikia1227731996.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/121k5z1227731996.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/135dks1227731996.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/14n7r81227731996.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/15803i1227731996.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/16nyj51227731996.tab")
+ }
>
> system("convert tmp/11iwd1227731996.ps tmp/11iwd1227731996.png")
> system("convert tmp/2l53j1227731996.ps tmp/2l53j1227731996.png")
> system("convert tmp/3rkht1227731996.ps tmp/3rkht1227731996.png")
> system("convert tmp/4yowk1227731996.ps tmp/4yowk1227731996.png")
> system("convert tmp/581i71227731996.ps tmp/581i71227731996.png")
> system("convert tmp/6dkd71227731996.ps tmp/6dkd71227731996.png")
> system("convert tmp/76lz21227731996.ps tmp/76lz21227731996.png")
> system("convert tmp/85q0q1227731996.ps tmp/85q0q1227731996.png")
> system("convert tmp/9qm981227731996.ps tmp/9qm981227731996.png")
> system("convert tmp/10nc731227731996.ps tmp/10nc731227731996.png")
>
>
> proc.time()
user system elapsed
4.022 2.503 4.408