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(99.2
+ ,11554.5
+ ,93.6
+ ,13182.1
+ ,104.2
+ ,14800.1
+ ,95.3
+ ,12150.7
+ ,102.7
+ ,14478.2
+ ,103.1
+ ,13253.9
+ ,100
+ ,12036.8
+ ,107.2
+ ,12653.2
+ ,107
+ ,14035.4
+ ,119
+ ,14571.4
+ ,110.4
+ ,15400.9
+ ,101.7
+ ,14283.2
+ ,102.4
+ ,14485.3
+ ,98.8
+ ,14196.3
+ ,105.6
+ ,15559.1
+ ,104.4
+ ,13767.4
+ ,106.3
+ ,14634
+ ,107.2
+ ,14381.1
+ ,108.5
+ ,12509.9
+ ,106.9
+ ,12122.3
+ ,114.2
+ ,13122.3
+ ,125.9
+ ,13908.7
+ ,110.6
+ ,13456.5
+ ,110.5
+ ,12441.6
+ ,106.7
+ ,12953
+ ,104.7
+ ,13057.2
+ ,107.4
+ ,14350.1
+ ,109.8
+ ,13830.2
+ ,103.4
+ ,13755.5
+ ,114.8
+ ,13574.4
+ ,114.3
+ ,12802.6
+ ,109.6
+ ,11737.3
+ ,118.3
+ ,13850.2
+ ,127.3
+ ,15081.8
+ ,112.3
+ ,13653.3
+ ,114.9
+ ,14019.1
+ ,108.2
+ ,13962
+ ,105.4
+ ,13768.7
+ ,122.1
+ ,14747.1
+ ,113.5
+ ,13858.1
+ ,110
+ ,13188
+ ,125.3
+ ,13693.1
+ ,114.3
+ ,12970
+ ,115.6
+ ,11392.8
+ ,127.1
+ ,13985.2
+ ,123
+ ,14994.7
+ ,122.2
+ ,13584.7
+ ,126.4
+ ,14257.8
+ ,112.7
+ ,13553.4
+ ,105.8
+ ,14007.3
+ ,120.9
+ ,16535.8
+ ,116.3
+ ,14721.4
+ ,115.7
+ ,13664.6
+ ,127.9
+ ,16405.9
+ ,108.3
+ ,13829.4
+ ,121.1
+ ,13735.6
+ ,128.6
+ ,15870.5
+ ,123.1
+ ,15962.4
+ ,127.7
+ ,15744.1
+ ,126.6
+ ,16083.7
+ ,118.4
+ ,14863.9
+ ,110
+ ,15533.1
+ ,129.6
+ ,17473.1
+ ,115.8
+ ,15925.5
+ ,125.9
+ ,15573.7
+ ,128.4
+ ,17495
+ ,114
+ ,14155.8
+ ,125.6
+ ,14913.9
+ ,128.5
+ ,17250.4
+ ,136.6
+ ,15879.8
+ ,133.1
+ ,17647.8
+ ,124.6
+ ,17749.9)
+ ,dim=c(2
+ ,72)
+ ,dimnames=list(c('Voeding'
+ ,'Invoer')
+ ,1:72))
> y <- array(NA,dim=c(2,72),dimnames=list(c('Voeding','Invoer'),1:72))
> 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
Voeding Invoer
1 99.2 11554.5
2 93.6 13182.1
3 104.2 14800.1
4 95.3 12150.7
5 102.7 14478.2
6 103.1 13253.9
7 100.0 12036.8
8 107.2 12653.2
9 107.0 14035.4
10 119.0 14571.4
11 110.4 15400.9
12 101.7 14283.2
13 102.4 14485.3
14 98.8 14196.3
15 105.6 15559.1
16 104.4 13767.4
17 106.3 14634.0
18 107.2 14381.1
19 108.5 12509.9
20 106.9 12122.3
21 114.2 13122.3
22 125.9 13908.7
23 110.6 13456.5
24 110.5 12441.6
25 106.7 12953.0
26 104.7 13057.2
27 107.4 14350.1
28 109.8 13830.2
29 103.4 13755.5
30 114.8 13574.4
31 114.3 12802.6
32 109.6 11737.3
33 118.3 13850.2
34 127.3 15081.8
35 112.3 13653.3
36 114.9 14019.1
37 108.2 13962.0
38 105.4 13768.7
39 122.1 14747.1
40 113.5 13858.1
41 110.0 13188.0
42 125.3 13693.1
43 114.3 12970.0
44 115.6 11392.8
45 127.1 13985.2
46 123.0 14994.7
47 122.2 13584.7
48 126.4 14257.8
49 112.7 13553.4
50 105.8 14007.3
51 120.9 16535.8
52 116.3 14721.4
53 115.7 13664.6
54 127.9 16405.9
55 108.3 13829.4
56 121.1 13735.6
57 128.6 15870.5
58 123.1 15962.4
59 127.7 15744.1
60 126.6 16083.7
61 118.4 14863.9
62 110.0 15533.1
63 129.6 17473.1
64 115.8 15925.5
65 125.9 15573.7
66 128.4 17495.0
67 114.0 14155.8
68 125.6 14913.9
69 128.5 17250.4
70 136.6 15879.8
71 133.1 17647.8
72 124.6 17749.9
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Invoer
53.328201 0.004265
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-15.9469 -6.3017 0.9068 5.6739 15.5480
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.333e+01 9.143e+00 5.833 1.54e-07 ***
Invoer 4.265e-03 6.367e-04 6.698 4.42e-09 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 7.815 on 70 degrees of freedom
Multiple R-squared: 0.3906, Adjusted R-squared: 0.3819
F-statistic: 44.86 on 1 and 70 DF, p-value: 4.415e-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,] 0.17280007 0.34560014 0.82719993
[2,] 0.11429009 0.22858018 0.88570991
[3,] 0.06186319 0.12372638 0.93813681
[4,] 0.11075232 0.22150463 0.88924768
[5,] 0.08615474 0.17230947 0.91384526
[6,] 0.30820857 0.61641715 0.69179143
[7,] 0.22995602 0.45991203 0.77004398
[8,] 0.21717807 0.43435614 0.78282193
[9,] 0.20524356 0.41048713 0.79475644
[10,] 0.25666151 0.51332302 0.74333849
[11,] 0.25966644 0.51933287 0.74033356
[12,] 0.22130977 0.44261954 0.77869023
[13,] 0.19924091 0.39848182 0.80075909
[14,] 0.17840639 0.35681277 0.82159361
[15,] 0.21029001 0.42058001 0.78970999
[16,] 0.20504520 0.41009040 0.79495480
[17,] 0.30423381 0.60846763 0.69576619
[18,] 0.78357042 0.43285915 0.21642958
[19,] 0.75366527 0.49266946 0.24633473
[20,] 0.72877369 0.54245262 0.27122631
[21,] 0.68329090 0.63341819 0.31670910
[22,] 0.65503952 0.68992095 0.34496048
[23,] 0.65412052 0.69175897 0.34587948
[24,] 0.62422198 0.75155603 0.37577802
[25,] 0.67902716 0.64194569 0.32097284
[26,] 0.68235323 0.63529353 0.31764677
[27,] 0.68540594 0.62918812 0.31459406
[28,] 0.64088483 0.71823034 0.35911517
[29,] 0.67442275 0.65115450 0.32557725
[30,] 0.82691369 0.34617262 0.17308631
[31,] 0.79623839 0.40752321 0.20376161
[32,] 0.76701413 0.46597174 0.23298587
[33,] 0.76838498 0.46323004 0.23161502
[34,] 0.81436486 0.37127028 0.18563514
[35,] 0.82476740 0.35046520 0.17523260
[36,] 0.79810525 0.40378950 0.20189475
[37,] 0.77745784 0.44508431 0.22254216
[38,] 0.86370992 0.27258015 0.13629008
[39,] 0.83553971 0.32892057 0.16446029
[40,] 0.84301010 0.31397981 0.15698990
[41,] 0.91892856 0.16214288 0.08107144
[42,] 0.90979504 0.18040991 0.09020496
[43,] 0.92556017 0.14887966 0.07443983
[44,] 0.95763035 0.08473930 0.04236965
[45,] 0.93737575 0.12524850 0.06262425
[46,] 0.95450986 0.09098029 0.04549014
[47,] 0.94423775 0.11152450 0.05576225
[48,] 0.92299103 0.15401793 0.07700897
[49,] 0.89131357 0.21737286 0.10868643
[50,] 0.86465574 0.27068853 0.13534426
[51,] 0.87921670 0.24156659 0.12078330
[52,] 0.85700123 0.28599754 0.14299877
[53,] 0.84077866 0.31844268 0.15922134
[54,] 0.78011363 0.43977274 0.21988637
[55,] 0.74608178 0.50783644 0.25391822
[56,] 0.67680491 0.64639017 0.32319509
[57,] 0.58071520 0.83856960 0.41928480
[58,] 0.75470971 0.49058057 0.24529029
[59,] 0.65450217 0.69099565 0.34549783
[60,] 0.72714086 0.54571829 0.27285914
[61,] 0.60907083 0.78185834 0.39092917
[62,] 0.46784545 0.93569089 0.53215455
[63,] 0.63379352 0.73241297 0.36620648
> postscript(file="/var/www/html/rcomp/tmp/1axb21229760739.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/2xrdd1229760739.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/35vzv1229760739.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/4kkxn1229760739.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/5ynnz1229760739.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 = 72
Frequency = 1
1 2 3 4 5 6
-3.40558190 -15.94693510 -12.24734644 -9.84824274 -12.37451430 -6.75314619
7 8 9 10 11 12
-4.66248448 -0.09129385 -6.18607045 3.52800836 -8.60962526 -12.54288252
13 14 15 16 17 18
-12.70479422 -15.07227329 -14.08431319 -7.64310986 -9.43896676 -7.46040432
19 20 21 22 23 24
1.81984889 1.87287697 4.90809864 13.25427696 -0.11719028 4.11113325
25 26 27 28 29 30
-1.86987439 -4.31426429 -7.12819619 -2.51093794 -8.59235900 3.57999236
31 32 33 34 35 36
6.37154827 6.21481663 5.90376649 9.65126550 0.74350135 1.78344543
37 38 39 40 41 42
-4.67303572 -6.64865407 5.87868681 1.07007475 0.42790270 13.57376317
43 44 45 46 47 48
5.65762438 13.68403276 14.12802142 5.72272770 10.93606514 12.26544285
49 50 51 52 53 54
1.56955270 -7.26623018 -2.94972219 0.18829161 4.09530935 4.60427252
55 56 57 58 59 60
-4.00752612 9.19251009 7.58763484 1.69570171 7.22670282 4.67838410
61 62 63 64 65 66
1.68056070 -9.57342896 1.75290109 -5.44692797 6.15342104 0.45950244
67 68 69 70 71 72
0.30045024 8.66732179 1.60266722 15.54797240 4.50784431 -4.42758956
> postscript(file="/var/www/html/rcomp/tmp/6kbf21229760739.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 = 72
Frequency = 1
lag(myerror, k = 1) myerror
0 -3.40558190 NA
1 -15.94693510 -3.40558190
2 -12.24734644 -15.94693510
3 -9.84824274 -12.24734644
4 -12.37451430 -9.84824274
5 -6.75314619 -12.37451430
6 -4.66248448 -6.75314619
7 -0.09129385 -4.66248448
8 -6.18607045 -0.09129385
9 3.52800836 -6.18607045
10 -8.60962526 3.52800836
11 -12.54288252 -8.60962526
12 -12.70479422 -12.54288252
13 -15.07227329 -12.70479422
14 -14.08431319 -15.07227329
15 -7.64310986 -14.08431319
16 -9.43896676 -7.64310986
17 -7.46040432 -9.43896676
18 1.81984889 -7.46040432
19 1.87287697 1.81984889
20 4.90809864 1.87287697
21 13.25427696 4.90809864
22 -0.11719028 13.25427696
23 4.11113325 -0.11719028
24 -1.86987439 4.11113325
25 -4.31426429 -1.86987439
26 -7.12819619 -4.31426429
27 -2.51093794 -7.12819619
28 -8.59235900 -2.51093794
29 3.57999236 -8.59235900
30 6.37154827 3.57999236
31 6.21481663 6.37154827
32 5.90376649 6.21481663
33 9.65126550 5.90376649
34 0.74350135 9.65126550
35 1.78344543 0.74350135
36 -4.67303572 1.78344543
37 -6.64865407 -4.67303572
38 5.87868681 -6.64865407
39 1.07007475 5.87868681
40 0.42790270 1.07007475
41 13.57376317 0.42790270
42 5.65762438 13.57376317
43 13.68403276 5.65762438
44 14.12802142 13.68403276
45 5.72272770 14.12802142
46 10.93606514 5.72272770
47 12.26544285 10.93606514
48 1.56955270 12.26544285
49 -7.26623018 1.56955270
50 -2.94972219 -7.26623018
51 0.18829161 -2.94972219
52 4.09530935 0.18829161
53 4.60427252 4.09530935
54 -4.00752612 4.60427252
55 9.19251009 -4.00752612
56 7.58763484 9.19251009
57 1.69570171 7.58763484
58 7.22670282 1.69570171
59 4.67838410 7.22670282
60 1.68056070 4.67838410
61 -9.57342896 1.68056070
62 1.75290109 -9.57342896
63 -5.44692797 1.75290109
64 6.15342104 -5.44692797
65 0.45950244 6.15342104
66 0.30045024 0.45950244
67 8.66732179 0.30045024
68 1.60266722 8.66732179
69 15.54797240 1.60266722
70 4.50784431 15.54797240
71 -4.42758956 4.50784431
72 NA -4.42758956
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -15.94693510 -3.40558190
[2,] -12.24734644 -15.94693510
[3,] -9.84824274 -12.24734644
[4,] -12.37451430 -9.84824274
[5,] -6.75314619 -12.37451430
[6,] -4.66248448 -6.75314619
[7,] -0.09129385 -4.66248448
[8,] -6.18607045 -0.09129385
[9,] 3.52800836 -6.18607045
[10,] -8.60962526 3.52800836
[11,] -12.54288252 -8.60962526
[12,] -12.70479422 -12.54288252
[13,] -15.07227329 -12.70479422
[14,] -14.08431319 -15.07227329
[15,] -7.64310986 -14.08431319
[16,] -9.43896676 -7.64310986
[17,] -7.46040432 -9.43896676
[18,] 1.81984889 -7.46040432
[19,] 1.87287697 1.81984889
[20,] 4.90809864 1.87287697
[21,] 13.25427696 4.90809864
[22,] -0.11719028 13.25427696
[23,] 4.11113325 -0.11719028
[24,] -1.86987439 4.11113325
[25,] -4.31426429 -1.86987439
[26,] -7.12819619 -4.31426429
[27,] -2.51093794 -7.12819619
[28,] -8.59235900 -2.51093794
[29,] 3.57999236 -8.59235900
[30,] 6.37154827 3.57999236
[31,] 6.21481663 6.37154827
[32,] 5.90376649 6.21481663
[33,] 9.65126550 5.90376649
[34,] 0.74350135 9.65126550
[35,] 1.78344543 0.74350135
[36,] -4.67303572 1.78344543
[37,] -6.64865407 -4.67303572
[38,] 5.87868681 -6.64865407
[39,] 1.07007475 5.87868681
[40,] 0.42790270 1.07007475
[41,] 13.57376317 0.42790270
[42,] 5.65762438 13.57376317
[43,] 13.68403276 5.65762438
[44,] 14.12802142 13.68403276
[45,] 5.72272770 14.12802142
[46,] 10.93606514 5.72272770
[47,] 12.26544285 10.93606514
[48,] 1.56955270 12.26544285
[49,] -7.26623018 1.56955270
[50,] -2.94972219 -7.26623018
[51,] 0.18829161 -2.94972219
[52,] 4.09530935 0.18829161
[53,] 4.60427252 4.09530935
[54,] -4.00752612 4.60427252
[55,] 9.19251009 -4.00752612
[56,] 7.58763484 9.19251009
[57,] 1.69570171 7.58763484
[58,] 7.22670282 1.69570171
[59,] 4.67838410 7.22670282
[60,] 1.68056070 4.67838410
[61,] -9.57342896 1.68056070
[62,] 1.75290109 -9.57342896
[63,] -5.44692797 1.75290109
[64,] 6.15342104 -5.44692797
[65,] 0.45950244 6.15342104
[66,] 0.30045024 0.45950244
[67,] 8.66732179 0.30045024
[68,] 1.60266722 8.66732179
[69,] 15.54797240 1.60266722
[70,] 4.50784431 15.54797240
[71,] -4.42758956 4.50784431
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -15.94693510 -3.40558190
2 -12.24734644 -15.94693510
3 -9.84824274 -12.24734644
4 -12.37451430 -9.84824274
5 -6.75314619 -12.37451430
6 -4.66248448 -6.75314619
7 -0.09129385 -4.66248448
8 -6.18607045 -0.09129385
9 3.52800836 -6.18607045
10 -8.60962526 3.52800836
11 -12.54288252 -8.60962526
12 -12.70479422 -12.54288252
13 -15.07227329 -12.70479422
14 -14.08431319 -15.07227329
15 -7.64310986 -14.08431319
16 -9.43896676 -7.64310986
17 -7.46040432 -9.43896676
18 1.81984889 -7.46040432
19 1.87287697 1.81984889
20 4.90809864 1.87287697
21 13.25427696 4.90809864
22 -0.11719028 13.25427696
23 4.11113325 -0.11719028
24 -1.86987439 4.11113325
25 -4.31426429 -1.86987439
26 -7.12819619 -4.31426429
27 -2.51093794 -7.12819619
28 -8.59235900 -2.51093794
29 3.57999236 -8.59235900
30 6.37154827 3.57999236
31 6.21481663 6.37154827
32 5.90376649 6.21481663
33 9.65126550 5.90376649
34 0.74350135 9.65126550
35 1.78344543 0.74350135
36 -4.67303572 1.78344543
37 -6.64865407 -4.67303572
38 5.87868681 -6.64865407
39 1.07007475 5.87868681
40 0.42790270 1.07007475
41 13.57376317 0.42790270
42 5.65762438 13.57376317
43 13.68403276 5.65762438
44 14.12802142 13.68403276
45 5.72272770 14.12802142
46 10.93606514 5.72272770
47 12.26544285 10.93606514
48 1.56955270 12.26544285
49 -7.26623018 1.56955270
50 -2.94972219 -7.26623018
51 0.18829161 -2.94972219
52 4.09530935 0.18829161
53 4.60427252 4.09530935
54 -4.00752612 4.60427252
55 9.19251009 -4.00752612
56 7.58763484 9.19251009
57 1.69570171 7.58763484
58 7.22670282 1.69570171
59 4.67838410 7.22670282
60 1.68056070 4.67838410
61 -9.57342896 1.68056070
62 1.75290109 -9.57342896
63 -5.44692797 1.75290109
64 6.15342104 -5.44692797
65 0.45950244 6.15342104
66 0.30045024 0.45950244
67 8.66732179 0.30045024
68 1.60266722 8.66732179
69 15.54797240 1.60266722
70 4.50784431 15.54797240
71 -4.42758956 4.50784431
> 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/7zmvh1229760739.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/8njnk1229760739.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/9ajvj1229760739.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/10nldy1229760739.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/11o46i1229760739.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/12bh4f1229760739.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/13xtet1229760740.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/1431m21229760740.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/15ug171229760740.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/16svpo1229760740.tab")
+ }
>
> system("convert tmp/1axb21229760739.ps tmp/1axb21229760739.png")
> system("convert tmp/2xrdd1229760739.ps tmp/2xrdd1229760739.png")
> system("convert tmp/35vzv1229760739.ps tmp/35vzv1229760739.png")
> system("convert tmp/4kkxn1229760739.ps tmp/4kkxn1229760739.png")
> system("convert tmp/5ynnz1229760739.ps tmp/5ynnz1229760739.png")
> system("convert tmp/6kbf21229760739.ps tmp/6kbf21229760739.png")
> system("convert tmp/7zmvh1229760739.ps tmp/7zmvh1229760739.png")
> system("convert tmp/8njnk1229760739.ps tmp/8njnk1229760739.png")
> system("convert tmp/9ajvj1229760739.ps tmp/9ajvj1229760739.png")
> system("convert tmp/10nldy1229760739.ps tmp/10nldy1229760739.png")
>
>
> proc.time()
user system elapsed
2.568 1.604 3.697