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(111.8
+ ,142
+ ,129.5
+ ,100.9
+ ,109
+ ,103.7
+ ,102.1
+ ,120.7
+ ,114.2
+ ,118
+ ,159.6
+ ,123.5
+ ,110.6
+ ,142.4
+ ,111.3
+ ,115.1
+ ,145.1
+ ,126.5
+ ,114.8
+ ,148.9
+ ,119.6
+ ,110.1
+ ,136.9
+ ,116.5
+ ,110.8
+ ,119.9
+ ,116.7
+ ,95.6
+ ,133.9
+ ,119.4
+ ,108.1
+ ,131
+ ,124
+ ,116
+ ,133.2
+ ,130.6
+ ,111.2
+ ,135
+ ,120.1
+ ,98.2
+ ,99.1
+ ,113.2
+ ,97.6
+ ,110.8
+ ,111.1
+ ,113.3
+ ,152.3
+ ,126
+ ,107
+ ,131.9
+ ,115.8
+ ,107.9
+ ,127.9
+ ,111
+ ,117.5
+ ,142
+ ,128.7
+ ,105.4
+ ,118.7
+ ,112.6
+ ,104.2
+ ,116.3
+ ,114.7
+ ,98
+ ,125.7
+ ,118.5
+ ,106.7
+ ,122.7
+ ,124.8
+ ,113.4
+ ,125.3
+ ,128.6
+ ,111.7
+ ,123.2
+ ,127
+ ,94.2
+ ,88.8
+ ,111.8
+ ,92.5
+ ,94.9
+ ,100.6
+ ,109.8
+ ,136.8
+ ,122.9
+ ,105.1
+ ,128.7
+ ,117.8
+ ,104.4
+ ,110.8
+ ,108.1
+ ,111.1
+ ,132.8
+ ,129.6
+ ,98.7
+ ,112
+ ,111.4
+ ,100.5
+ ,104.5
+ ,110
+ ,93.7
+ ,112
+ ,115.2
+ ,103.2
+ ,110.6
+ ,118.8
+ ,104.1
+ ,107.2
+ ,116.2
+ ,106.9
+ ,116.2
+ ,126.3
+ ,89.2
+ ,85.7
+ ,106.7
+ ,88.7
+ ,94.2
+ ,96.5
+ ,110.7
+ ,127.2
+ ,119.1
+ ,98.8
+ ,108.9
+ ,109.6
+ ,102.5
+ ,111.9
+ ,110.3
+ ,101.8
+ ,126.3
+ ,118.8
+ ,96
+ ,105.9
+ ,104.5
+ ,98.3
+ ,101.3
+ ,107.7
+ ,94
+ ,105.5
+ ,127.7
+ ,105.1
+ ,106.3
+ ,118.5
+ ,114
+ ,117.3
+ ,120.1
+ ,115.5
+ ,110.9
+ ,127.4
+ ,94.3
+ ,85.4
+ ,107.8
+ ,100.8
+ ,81.9
+ ,106.5
+ ,111.2
+ ,121.5
+ ,124.6
+ ,103.4
+ ,106.3
+ ,101.9
+ ,106.7
+ ,111.8
+ ,106.5
+ ,112.2
+ ,122.8
+ ,119.4
+ ,100.7
+ ,101.8
+ ,103.3
+ ,99
+ ,92.2
+ ,99.6
+ ,91.5
+ ,106.3
+ ,120.9
+ ,102.7
+ ,103
+ ,111.7
+ ,111.4
+ ,97.7
+ ,123.9)
+ ,dim=c(3
+ ,60)
+ ,dimnames=list(c('Interm.'
+ ,'Invest.'
+ ,'Cons.')
+ ,1:60))
> y <- array(NA,dim=c(3,60),dimnames=list(c('Interm.','Invest.','Cons.'),1:60))
> 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 = '3'
> #'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
Cons. Interm. Invest.
1 129.5 111.8 142.0
2 103.7 100.9 109.0
3 114.2 102.1 120.7
4 123.5 118.0 159.6
5 111.3 110.6 142.4
6 126.5 115.1 145.1
7 119.6 114.8 148.9
8 116.5 110.1 136.9
9 116.7 110.8 119.9
10 119.4 95.6 133.9
11 124.0 108.1 131.0
12 130.6 116.0 133.2
13 120.1 111.2 135.0
14 113.2 98.2 99.1
15 111.1 97.6 110.8
16 126.0 113.3 152.3
17 115.8 107.0 131.9
18 111.0 107.9 127.9
19 128.7 117.5 142.0
20 112.6 105.4 118.7
21 114.7 104.2 116.3
22 118.5 98.0 125.7
23 124.8 106.7 122.7
24 128.6 113.4 125.3
25 127.0 111.7 123.2
26 111.8 94.2 88.8
27 100.6 92.5 94.9
28 122.9 109.8 136.8
29 117.8 105.1 128.7
30 108.1 104.4 110.8
31 129.6 111.1 132.8
32 111.4 98.7 112.0
33 110.0 100.5 104.5
34 115.2 93.7 112.0
35 118.8 103.2 110.6
36 116.2 104.1 107.2
37 126.3 106.9 116.2
38 106.7 89.2 85.7
39 96.5 88.7 94.2
40 119.1 110.7 127.2
41 109.6 98.8 108.9
42 110.3 102.5 111.9
43 118.8 101.8 126.3
44 104.5 96.0 105.9
45 107.7 98.3 101.3
46 127.7 94.0 105.5
47 118.5 105.1 106.3
48 120.1 114.0 117.3
49 127.4 115.5 110.9
50 107.8 94.3 85.4
51 106.5 100.8 81.9
52 124.6 111.2 121.5
53 101.9 103.4 106.3
54 106.5 106.7 111.8
55 119.4 112.2 122.8
56 103.3 100.7 101.8
57 99.6 99.0 92.2
58 120.9 91.5 106.3
59 111.7 102.7 103.0
60 123.9 111.4 97.7
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Interm. Invest.
47.4205 0.4751 0.1620
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-11.8731 -3.9979 -0.5693 4.4833 18.5226
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 47.42046 12.10942 3.916 0.000243 ***
Interm. 0.47512 0.15870 2.994 0.004070 **
Invest. 0.16205 0.06852 2.365 0.021461 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 6.311 on 57 degrees of freedom
Multiple R-squared: 0.4815, Adjusted R-squared: 0.4633
F-statistic: 26.47 on 2 and 57 DF, p-value: 7.406e-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.837407304 0.32518539 0.1625927
[2,] 0.766938386 0.46612323 0.2330616
[3,] 0.664873110 0.67025378 0.3351269
[4,] 0.545874389 0.90825122 0.4541256
[5,] 0.539635280 0.92072944 0.4603647
[6,] 0.549944838 0.90011032 0.4500552
[7,] 0.667148077 0.66570385 0.3328519
[8,] 0.574638040 0.85072392 0.4253620
[9,] 0.505470494 0.98905901 0.4945295
[10,] 0.410314954 0.82062991 0.5896850
[11,] 0.331995144 0.66399029 0.6680049
[12,] 0.276795284 0.55359057 0.7232047
[13,] 0.325314710 0.65062942 0.6746853
[14,] 0.283695459 0.56739092 0.7163045
[15,] 0.239170874 0.47834175 0.7608291
[16,] 0.180954158 0.36190832 0.8190458
[17,] 0.154483525 0.30896705 0.8455165
[18,] 0.176079194 0.35215839 0.8239208
[19,] 0.195674774 0.39134955 0.8043252
[20,] 0.191743924 0.38348785 0.8082561
[21,] 0.167406639 0.33481328 0.8325934
[22,] 0.172775756 0.34555151 0.8272242
[23,] 0.130944601 0.26188920 0.8690554
[24,] 0.096968707 0.19393741 0.9030313
[25,] 0.110251428 0.22050286 0.8897486
[26,] 0.118838195 0.23767639 0.8811618
[27,] 0.086209502 0.17241900 0.9137905
[28,] 0.062286536 0.12457307 0.9377135
[29,] 0.053646336 0.10729267 0.9463537
[30,] 0.041727701 0.08345540 0.9582723
[31,] 0.027647072 0.05529414 0.9723529
[32,] 0.038783835 0.07756767 0.9612162
[33,] 0.027773255 0.05554651 0.9722267
[34,] 0.035960468 0.07192094 0.9640395
[35,] 0.024000913 0.04800183 0.9759991
[36,] 0.016117679 0.03223536 0.9838823
[37,] 0.012289275 0.02457855 0.9877107
[38,] 0.007618015 0.01523603 0.9923820
[39,] 0.007878819 0.01575764 0.9921212
[40,] 0.005268799 0.01053760 0.9947312
[41,] 0.082354172 0.16470834 0.9176458
[42,] 0.061004764 0.12200953 0.9389952
[43,] 0.037271198 0.07454240 0.9627288
[44,] 0.041250533 0.08250107 0.9587495
[45,] 0.025913476 0.05182695 0.9740865
[46,] 0.015623921 0.03124784 0.9843761
[47,] 0.010981527 0.02196305 0.9890185
[48,] 0.020475462 0.04095092 0.9795245
[49,] 0.027306016 0.05461203 0.9726940
> postscript(file="/var/www/html/rcomp/tmp/1cch21229415703.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/28jcn1229415703.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/3rz3a1229415703.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/4yi241229415703.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/5l0ek1229415703.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 = 60
Frequency = 1
1 2 3 4 5 6
5.95083839 -9.32287220 -1.28894512 -5.84689516 -11.74383741 0.88060638
7 8 9 10 11 12
-6.49263036 -5.41502868 -2.79284079 4.86032729 3.99127610 6.48133927
13 14 15 16 17 18
-2.02977300 3.06419701 -0.64666240 0.06909338 -3.83193423 -8.41135959
19 20 21 22 23 24
2.44266228 -4.13274585 -1.07369468 4.14881447 6.80141856 6.99680592
25 26 27 28 29 30
6.54480280 5.23373856 -6.14703642 1.14371143 -0.41066375 -6.87746898
31 32 33 34 35 36
7.87423863 -1.06374729 -2.10362069 5.11184579 4.42508243 1.94842986
37 38 39 40 41 42
9.25968960 3.01167222 -8.32815393 -1.52825998 -2.40891857 -3.95299350
43 44 45 46 47 48
2.54613652 -5.69245040 -2.83981462 18.52260497 3.91915205 -0.49190246
49 50 51 52 53 54
7.13250985 1.73718088 -2.08393140 4.65783920 -11.87314630 -9.73228715
55 56 57 58 59 60
-1.22793837 -8.46112197 -9.79778498 12.78076523 -1.20581362 7.71949477
> postscript(file="/var/www/html/rcomp/tmp/6e2gd1229415703.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 = 60
Frequency = 1
lag(myerror, k = 1) myerror
0 5.95083839 NA
1 -9.32287220 5.95083839
2 -1.28894512 -9.32287220
3 -5.84689516 -1.28894512
4 -11.74383741 -5.84689516
5 0.88060638 -11.74383741
6 -6.49263036 0.88060638
7 -5.41502868 -6.49263036
8 -2.79284079 -5.41502868
9 4.86032729 -2.79284079
10 3.99127610 4.86032729
11 6.48133927 3.99127610
12 -2.02977300 6.48133927
13 3.06419701 -2.02977300
14 -0.64666240 3.06419701
15 0.06909338 -0.64666240
16 -3.83193423 0.06909338
17 -8.41135959 -3.83193423
18 2.44266228 -8.41135959
19 -4.13274585 2.44266228
20 -1.07369468 -4.13274585
21 4.14881447 -1.07369468
22 6.80141856 4.14881447
23 6.99680592 6.80141856
24 6.54480280 6.99680592
25 5.23373856 6.54480280
26 -6.14703642 5.23373856
27 1.14371143 -6.14703642
28 -0.41066375 1.14371143
29 -6.87746898 -0.41066375
30 7.87423863 -6.87746898
31 -1.06374729 7.87423863
32 -2.10362069 -1.06374729
33 5.11184579 -2.10362069
34 4.42508243 5.11184579
35 1.94842986 4.42508243
36 9.25968960 1.94842986
37 3.01167222 9.25968960
38 -8.32815393 3.01167222
39 -1.52825998 -8.32815393
40 -2.40891857 -1.52825998
41 -3.95299350 -2.40891857
42 2.54613652 -3.95299350
43 -5.69245040 2.54613652
44 -2.83981462 -5.69245040
45 18.52260497 -2.83981462
46 3.91915205 18.52260497
47 -0.49190246 3.91915205
48 7.13250985 -0.49190246
49 1.73718088 7.13250985
50 -2.08393140 1.73718088
51 4.65783920 -2.08393140
52 -11.87314630 4.65783920
53 -9.73228715 -11.87314630
54 -1.22793837 -9.73228715
55 -8.46112197 -1.22793837
56 -9.79778498 -8.46112197
57 12.78076523 -9.79778498
58 -1.20581362 12.78076523
59 7.71949477 -1.20581362
60 NA 7.71949477
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -9.32287220 5.95083839
[2,] -1.28894512 -9.32287220
[3,] -5.84689516 -1.28894512
[4,] -11.74383741 -5.84689516
[5,] 0.88060638 -11.74383741
[6,] -6.49263036 0.88060638
[7,] -5.41502868 -6.49263036
[8,] -2.79284079 -5.41502868
[9,] 4.86032729 -2.79284079
[10,] 3.99127610 4.86032729
[11,] 6.48133927 3.99127610
[12,] -2.02977300 6.48133927
[13,] 3.06419701 -2.02977300
[14,] -0.64666240 3.06419701
[15,] 0.06909338 -0.64666240
[16,] -3.83193423 0.06909338
[17,] -8.41135959 -3.83193423
[18,] 2.44266228 -8.41135959
[19,] -4.13274585 2.44266228
[20,] -1.07369468 -4.13274585
[21,] 4.14881447 -1.07369468
[22,] 6.80141856 4.14881447
[23,] 6.99680592 6.80141856
[24,] 6.54480280 6.99680592
[25,] 5.23373856 6.54480280
[26,] -6.14703642 5.23373856
[27,] 1.14371143 -6.14703642
[28,] -0.41066375 1.14371143
[29,] -6.87746898 -0.41066375
[30,] 7.87423863 -6.87746898
[31,] -1.06374729 7.87423863
[32,] -2.10362069 -1.06374729
[33,] 5.11184579 -2.10362069
[34,] 4.42508243 5.11184579
[35,] 1.94842986 4.42508243
[36,] 9.25968960 1.94842986
[37,] 3.01167222 9.25968960
[38,] -8.32815393 3.01167222
[39,] -1.52825998 -8.32815393
[40,] -2.40891857 -1.52825998
[41,] -3.95299350 -2.40891857
[42,] 2.54613652 -3.95299350
[43,] -5.69245040 2.54613652
[44,] -2.83981462 -5.69245040
[45,] 18.52260497 -2.83981462
[46,] 3.91915205 18.52260497
[47,] -0.49190246 3.91915205
[48,] 7.13250985 -0.49190246
[49,] 1.73718088 7.13250985
[50,] -2.08393140 1.73718088
[51,] 4.65783920 -2.08393140
[52,] -11.87314630 4.65783920
[53,] -9.73228715 -11.87314630
[54,] -1.22793837 -9.73228715
[55,] -8.46112197 -1.22793837
[56,] -9.79778498 -8.46112197
[57,] 12.78076523 -9.79778498
[58,] -1.20581362 12.78076523
[59,] 7.71949477 -1.20581362
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -9.32287220 5.95083839
2 -1.28894512 -9.32287220
3 -5.84689516 -1.28894512
4 -11.74383741 -5.84689516
5 0.88060638 -11.74383741
6 -6.49263036 0.88060638
7 -5.41502868 -6.49263036
8 -2.79284079 -5.41502868
9 4.86032729 -2.79284079
10 3.99127610 4.86032729
11 6.48133927 3.99127610
12 -2.02977300 6.48133927
13 3.06419701 -2.02977300
14 -0.64666240 3.06419701
15 0.06909338 -0.64666240
16 -3.83193423 0.06909338
17 -8.41135959 -3.83193423
18 2.44266228 -8.41135959
19 -4.13274585 2.44266228
20 -1.07369468 -4.13274585
21 4.14881447 -1.07369468
22 6.80141856 4.14881447
23 6.99680592 6.80141856
24 6.54480280 6.99680592
25 5.23373856 6.54480280
26 -6.14703642 5.23373856
27 1.14371143 -6.14703642
28 -0.41066375 1.14371143
29 -6.87746898 -0.41066375
30 7.87423863 -6.87746898
31 -1.06374729 7.87423863
32 -2.10362069 -1.06374729
33 5.11184579 -2.10362069
34 4.42508243 5.11184579
35 1.94842986 4.42508243
36 9.25968960 1.94842986
37 3.01167222 9.25968960
38 -8.32815393 3.01167222
39 -1.52825998 -8.32815393
40 -2.40891857 -1.52825998
41 -3.95299350 -2.40891857
42 2.54613652 -3.95299350
43 -5.69245040 2.54613652
44 -2.83981462 -5.69245040
45 18.52260497 -2.83981462
46 3.91915205 18.52260497
47 -0.49190246 3.91915205
48 7.13250985 -0.49190246
49 1.73718088 7.13250985
50 -2.08393140 1.73718088
51 4.65783920 -2.08393140
52 -11.87314630 4.65783920
53 -9.73228715 -11.87314630
54 -1.22793837 -9.73228715
55 -8.46112197 -1.22793837
56 -9.79778498 -8.46112197
57 12.78076523 -9.79778498
58 -1.20581362 12.78076523
59 7.71949477 -1.20581362
> 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/7t77n1229415703.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/8ja361229415703.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/9ygv41229415703.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/10kmm31229415703.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/11z1lu1229415703.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/120azi1229415704.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/13d8wk1229415704.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/147rq21229415704.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/15o9le1229415704.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/16uyjq1229415704.tab")
+ }
>
> system("convert tmp/1cch21229415703.ps tmp/1cch21229415703.png")
> system("convert tmp/28jcn1229415703.ps tmp/28jcn1229415703.png")
> system("convert tmp/3rz3a1229415703.ps tmp/3rz3a1229415703.png")
> system("convert tmp/4yi241229415703.ps tmp/4yi241229415703.png")
> system("convert tmp/5l0ek1229415703.ps tmp/5l0ek1229415703.png")
> system("convert tmp/6e2gd1229415703.ps tmp/6e2gd1229415703.png")
> system("convert tmp/7t77n1229415703.ps tmp/7t77n1229415703.png")
> system("convert tmp/8ja361229415703.ps tmp/8ja361229415703.png")
> system("convert tmp/9ygv41229415703.ps tmp/9ygv41229415703.png")
> system("convert tmp/10kmm31229415703.ps tmp/10kmm31229415703.png")
>
>
> proc.time()
user system elapsed
2.554 1.611 3.092