R version 2.9.0 (2009-04-17)
Copyright (C) 2009 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(114.08,136.49,112.95,142.62,135.31,141.71,134.31,149.51,133.03,147.39,140.11,131.96,124.69,136.38,131.68,127.34,150.95,133.85,137.26,125.14,130.51,141.25,143.15,149.32,118.01,120.92,122.56,134.85,147.97,131.93,135.74,134.22,151.62,143.07,154.82,145.37,145.59,134.32,147.12,126.31,175.86,162.21,140.66,124.09,152.69,153.91,154.38,154.34,132.45,138.70,136.44,150.98,153.24,146.39,154.11,178.30,155.93,168.23,142.53,162.52,148.73,158.86,147.73,152.17,166.79,171.01,144.30,171.49,156.07,189.62,161.70,177.46,152.10,179.98,140.45,156.96,155.56,167.89,174.53,194.78,167.16,192.78,159.48,165.06,173.22,196.60,176.13,151.64,180.31,187.02,185.84,210.99,169.43,219.08,195.25,235.68,174.99,241.44,156.42,187.46,182.08,229.57,182.00,208.44,153.28,215.09,136.72,217.00,130.19,171.08,132.04,178.41,143.89,196.34,133.38,172.11,127.98,154.93,150.45,182.26,133.55,181.74),dim=c(2,61),dimnames=list(c('InvoerEU','InvoerAM'),1:61))
> y <- array(NA,dim=c(2,61),dimnames=list(c('InvoerEU','InvoerAM'),1:61))
> for (i in 1:dim(x)[1])
+ {
+ for (j in 1:dim(x)[2])
+ {
+ y[i,j] <- as.numeric(x[i,j])
+ }
+ }
> par3 = 'Linear Trend'
> par2 = 'Include Monthly Dummies'
> par1 = '1'
> #'GNU S' R Code compiled by R2WASP v. 1.0.44 ()
> #Author: Prof. Dr. P. Wessa
> #To cite this work: AUTHOR(S), (YEAR), YOUR SOFTWARE TITLE (vNUMBER) in Free Statistics Software (v$_version), Office for Research Development and Education, URL http://www.wessa.net/rwasp_YOURPAGE.wasp/
> #Source of accompanying publication: Office for Research, Development, and Education
> #Technical description: Write here your technical program description (don't use hard returns!)
> library(lattice)
> library(lmtest)
Loading required package: zoo
Attaching package: 'zoo'
The following object(s) are masked from package:base :
as.Date.numeric
> n25 <- 25 #minimum number of obs. for Goldfeld-Quandt test
> par1 <- as.numeric(par1)
> x <- t(y)
> k <- length(x[1,])
> n <- length(x[,1])
> x1 <- cbind(x[,par1], x[,1:k!=par1])
> mycolnames <- c(colnames(x)[par1], colnames(x)[1:k!=par1])
> colnames(x1) <- mycolnames #colnames(x)[par1]
> x <- x1
> if (par3 == 'First Differences'){
+ x2 <- array(0, dim=c(n-1,k), dimnames=list(1:(n-1), paste('(1-B)',colnames(x),sep='')))
+ for (i in 1:n-1) {
+ for (j in 1:k) {
+ x2[i,j] <- x[i+1,j] - x[i,j]
+ }
+ }
+ x <- x2
+ }
> if (par2 == 'Include Monthly Dummies'){
+ x2 <- array(0, dim=c(n,11), dimnames=list(1:n, paste('M', seq(1:11), sep ='')))
+ for (i in 1:11){
+ x2[seq(i,n,12),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> if (par2 == 'Include Quarterly Dummies'){
+ x2 <- array(0, dim=c(n,3), dimnames=list(1:n, paste('Q', seq(1:3), sep ='')))
+ for (i in 1:3){
+ x2[seq(i,n,4),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> k <- length(x[1,])
> if (par3 == 'Linear Trend'){
+ x <- cbind(x, c(1:n))
+ colnames(x)[k+1] <- 't'
+ }
> x
InvoerEU InvoerAM M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 114.08 136.49 1 0 0 0 0 0 0 0 0 0 0 1
2 112.95 142.62 0 1 0 0 0 0 0 0 0 0 0 2
3 135.31 141.71 0 0 1 0 0 0 0 0 0 0 0 3
4 134.31 149.51 0 0 0 1 0 0 0 0 0 0 0 4
5 133.03 147.39 0 0 0 0 1 0 0 0 0 0 0 5
6 140.11 131.96 0 0 0 0 0 1 0 0 0 0 0 6
7 124.69 136.38 0 0 0 0 0 0 1 0 0 0 0 7
8 131.68 127.34 0 0 0 0 0 0 0 1 0 0 0 8
9 150.95 133.85 0 0 0 0 0 0 0 0 1 0 0 9
10 137.26 125.14 0 0 0 0 0 0 0 0 0 1 0 10
11 130.51 141.25 0 0 0 0 0 0 0 0 0 0 1 11
12 143.15 149.32 0 0 0 0 0 0 0 0 0 0 0 12
13 118.01 120.92 1 0 0 0 0 0 0 0 0 0 0 13
14 122.56 134.85 0 1 0 0 0 0 0 0 0 0 0 14
15 147.97 131.93 0 0 1 0 0 0 0 0 0 0 0 15
16 135.74 134.22 0 0 0 1 0 0 0 0 0 0 0 16
17 151.62 143.07 0 0 0 0 1 0 0 0 0 0 0 17
18 154.82 145.37 0 0 0 0 0 1 0 0 0 0 0 18
19 145.59 134.32 0 0 0 0 0 0 1 0 0 0 0 19
20 147.12 126.31 0 0 0 0 0 0 0 1 0 0 0 20
21 175.86 162.21 0 0 0 0 0 0 0 0 1 0 0 21
22 140.66 124.09 0 0 0 0 0 0 0 0 0 1 0 22
23 152.69 153.91 0 0 0 0 0 0 0 0 0 0 1 23
24 154.38 154.34 0 0 0 0 0 0 0 0 0 0 0 24
25 132.45 138.70 1 0 0 0 0 0 0 0 0 0 0 25
26 136.44 150.98 0 1 0 0 0 0 0 0 0 0 0 26
27 153.24 146.39 0 0 1 0 0 0 0 0 0 0 0 27
28 154.11 178.30 0 0 0 1 0 0 0 0 0 0 0 28
29 155.93 168.23 0 0 0 0 1 0 0 0 0 0 0 29
30 142.53 162.52 0 0 0 0 0 1 0 0 0 0 0 30
31 148.73 158.86 0 0 0 0 0 0 1 0 0 0 0 31
32 147.73 152.17 0 0 0 0 0 0 0 1 0 0 0 32
33 166.79 171.01 0 0 0 0 0 0 0 0 1 0 0 33
34 144.30 171.49 0 0 0 0 0 0 0 0 0 1 0 34
35 156.07 189.62 0 0 0 0 0 0 0 0 0 0 1 35
36 161.70 177.46 0 0 0 0 0 0 0 0 0 0 0 36
37 152.10 179.98 1 0 0 0 0 0 0 0 0 0 0 37
38 140.45 156.96 0 1 0 0 0 0 0 0 0 0 0 38
39 155.56 167.89 0 0 1 0 0 0 0 0 0 0 0 39
40 174.53 194.78 0 0 0 1 0 0 0 0 0 0 0 40
41 167.16 192.78 0 0 0 0 1 0 0 0 0 0 0 41
42 159.48 165.06 0 0 0 0 0 1 0 0 0 0 0 42
43 173.22 196.60 0 0 0 0 0 0 1 0 0 0 0 43
44 176.13 151.64 0 0 0 0 0 0 0 1 0 0 0 44
45 180.31 187.02 0 0 0 0 0 0 0 0 1 0 0 45
46 185.84 210.99 0 0 0 0 0 0 0 0 0 1 0 46
47 169.43 219.08 0 0 0 0 0 0 0 0 0 0 1 47
48 195.25 235.68 0 0 0 0 0 0 0 0 0 0 0 48
49 174.99 241.44 1 0 0 0 0 0 0 0 0 0 0 49
50 156.42 187.46 0 1 0 0 0 0 0 0 0 0 0 50
51 182.08 229.57 0 0 1 0 0 0 0 0 0 0 0 51
52 182.00 208.44 0 0 0 1 0 0 0 0 0 0 0 52
53 153.28 215.09 0 0 0 0 1 0 0 0 0 0 0 53
54 136.72 217.00 0 0 0 0 0 1 0 0 0 0 0 54
55 130.19 171.08 0 0 0 0 0 0 1 0 0 0 0 55
56 132.04 178.41 0 0 0 0 0 0 0 1 0 0 0 56
57 143.89 196.34 0 0 0 0 0 0 0 0 1 0 0 57
58 133.38 172.11 0 0 0 0 0 0 0 0 0 1 0 58
59 127.98 154.93 0 0 0 0 0 0 0 0 0 0 1 59
60 150.45 182.26 0 0 0 0 0 0 0 0 0 0 0 60
61 133.55 181.74 1 0 0 0 0 0 0 0 0 0 0 61
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) InvoerAM M1 M2 M3 M4
82.6839 0.4673 -18.0509 -17.0174 0.0396 -2.9589
M5 M6 M7 M8 M9 M10
-6.8565 -7.9970 -7.7807 0.5693 6.6422 -4.1152
M11 t
-10.0455 -0.1588
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-30.787 -7.681 1.813 8.023 29.009
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 82.68385 14.78653 5.592 1.11e-06 ***
InvoerAM 0.46726 0.09868 4.735 2.05e-05 ***
M1 -18.05090 8.13176 -2.220 0.0313 *
M2 -17.01736 8.58962 -1.981 0.0534 .
M3 0.03960 8.50755 0.005 0.9963
M4 -2.95888 8.49490 -0.348 0.7292
M5 -6.85649 8.48233 -0.808 0.4230
M6 -7.99702 8.50286 -0.941 0.3518
M7 -7.78073 8.57245 -0.908 0.3687
M8 0.56927 8.87930 0.064 0.9492
M9 6.64216 8.47414 0.784 0.4371
M10 -4.11521 8.60581 -0.478 0.6347
M11 -10.04548 8.47652 -1.185 0.2419
t -0.15881 0.16343 -0.972 0.3361
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 13.36 on 47 degrees of freedom
Multiple R-squared: 0.5931, Adjusted R-squared: 0.4805
F-statistic: 5.269 on 13 and 47 DF, p-value: 1.078e-05
> 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,] 5.290334e-03 1.058067e-02 0.9947097
[2,] 1.706180e-02 3.412361e-02 0.9829382
[3,] 1.215686e-02 2.431371e-02 0.9878431
[4,] 3.598910e-03 7.197821e-03 0.9964011
[5,] 1.041638e-03 2.083277e-03 0.9989584
[6,] 1.181092e-03 2.362183e-03 0.9988189
[7,] 3.936977e-04 7.873954e-04 0.9996063
[8,] 1.444225e-04 2.888449e-04 0.9998556
[9,] 8.871371e-05 1.774274e-04 0.9999113
[10,] 3.695891e-05 7.391781e-05 0.9999630
[11,] 3.649655e-05 7.299310e-05 0.9999635
[12,] 4.598220e-05 9.196439e-05 0.9999540
[13,] 2.750243e-05 5.500485e-05 0.9999725
[14,] 1.031390e-03 2.062781e-03 0.9989686
[15,] 4.504514e-04 9.009029e-04 0.9995495
[16,] 3.555462e-04 7.110925e-04 0.9996445
[17,] 2.845641e-04 5.691282e-04 0.9997154
[18,] 4.102257e-04 8.204513e-04 0.9995898
[19,] 2.913001e-04 5.826001e-04 0.9997087
[20,] 3.329291e-04 6.658581e-04 0.9996671
[21,] 7.926089e-04 1.585218e-03 0.9992074
[22,] 1.409581e-03 2.819161e-03 0.9985904
[23,] 1.230606e-02 2.461211e-02 0.9876939
[24,] 4.034385e-01 8.068769e-01 0.5965615
[25,] 7.423143e-01 5.153713e-01 0.2576857
[26,] 8.956508e-01 2.086984e-01 0.1043492
[27,] 8.193443e-01 3.613114e-01 0.1806557
[28,] 7.411718e-01 5.176564e-01 0.2588282
> postscript(file="/var/www/html/rcomp/tmp/111u91259087225.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/2ro651259087225.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/3kksh1259087225.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/4qn0e1259087225.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/5cy6x1259087225.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 = 61
Frequency = 1
1 2 3 4 5 6
-14.1708420 -19.0398810 -13.1528185 -14.6401740 -10.8731540 4.7160542
7 8 9 10 11 12
-12.8267221 -9.8038500 0.5101946 1.8062297 -6.3822868 -7.3997630
13 14 15 16 17 18
-1.0597898 -3.8934786 5.9827821 -4.1599554 11.6411916 15.0658302
19 20 21 22 23 24
10.9416096 8.0232010 14.0743916 7.6026260 11.7879363 3.3903480
25 26 27 28 29 30
6.9780478 4.3553425 6.4019321 -4.4811297 6.1006296 -3.3319566
31 32 33 34 35 36
4.5207505 -1.5444450 2.7982492 -8.9998609 0.3877517 1.8130021
37 38 39 40 41 42
9.2452093 7.4768812 0.5815520 10.1441494 7.7650979 14.3369661
43 44 45 46 47 48
13.2820224 29.0089746 10.7431419 15.9890284 1.8879596 10.0647317
49 50 51 52 53 54
5.3230073 11.1011358 0.1865522 13.1371098 -14.6337651 -30.7868939
55 56 57 58 59 60
-15.9176604 -25.6838805 -28.1259772 -16.3980233 -7.6813609 -7.8683188
61
-6.3156327
> postscript(file="/var/www/html/rcomp/tmp/6bffd1259087225.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 = 61
Frequency = 1
lag(myerror, k = 1) myerror
0 -14.1708420 NA
1 -19.0398810 -14.1708420
2 -13.1528185 -19.0398810
3 -14.6401740 -13.1528185
4 -10.8731540 -14.6401740
5 4.7160542 -10.8731540
6 -12.8267221 4.7160542
7 -9.8038500 -12.8267221
8 0.5101946 -9.8038500
9 1.8062297 0.5101946
10 -6.3822868 1.8062297
11 -7.3997630 -6.3822868
12 -1.0597898 -7.3997630
13 -3.8934786 -1.0597898
14 5.9827821 -3.8934786
15 -4.1599554 5.9827821
16 11.6411916 -4.1599554
17 15.0658302 11.6411916
18 10.9416096 15.0658302
19 8.0232010 10.9416096
20 14.0743916 8.0232010
21 7.6026260 14.0743916
22 11.7879363 7.6026260
23 3.3903480 11.7879363
24 6.9780478 3.3903480
25 4.3553425 6.9780478
26 6.4019321 4.3553425
27 -4.4811297 6.4019321
28 6.1006296 -4.4811297
29 -3.3319566 6.1006296
30 4.5207505 -3.3319566
31 -1.5444450 4.5207505
32 2.7982492 -1.5444450
33 -8.9998609 2.7982492
34 0.3877517 -8.9998609
35 1.8130021 0.3877517
36 9.2452093 1.8130021
37 7.4768812 9.2452093
38 0.5815520 7.4768812
39 10.1441494 0.5815520
40 7.7650979 10.1441494
41 14.3369661 7.7650979
42 13.2820224 14.3369661
43 29.0089746 13.2820224
44 10.7431419 29.0089746
45 15.9890284 10.7431419
46 1.8879596 15.9890284
47 10.0647317 1.8879596
48 5.3230073 10.0647317
49 11.1011358 5.3230073
50 0.1865522 11.1011358
51 13.1371098 0.1865522
52 -14.6337651 13.1371098
53 -30.7868939 -14.6337651
54 -15.9176604 -30.7868939
55 -25.6838805 -15.9176604
56 -28.1259772 -25.6838805
57 -16.3980233 -28.1259772
58 -7.6813609 -16.3980233
59 -7.8683188 -7.6813609
60 -6.3156327 -7.8683188
61 NA -6.3156327
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -19.0398810 -14.1708420
[2,] -13.1528185 -19.0398810
[3,] -14.6401740 -13.1528185
[4,] -10.8731540 -14.6401740
[5,] 4.7160542 -10.8731540
[6,] -12.8267221 4.7160542
[7,] -9.8038500 -12.8267221
[8,] 0.5101946 -9.8038500
[9,] 1.8062297 0.5101946
[10,] -6.3822868 1.8062297
[11,] -7.3997630 -6.3822868
[12,] -1.0597898 -7.3997630
[13,] -3.8934786 -1.0597898
[14,] 5.9827821 -3.8934786
[15,] -4.1599554 5.9827821
[16,] 11.6411916 -4.1599554
[17,] 15.0658302 11.6411916
[18,] 10.9416096 15.0658302
[19,] 8.0232010 10.9416096
[20,] 14.0743916 8.0232010
[21,] 7.6026260 14.0743916
[22,] 11.7879363 7.6026260
[23,] 3.3903480 11.7879363
[24,] 6.9780478 3.3903480
[25,] 4.3553425 6.9780478
[26,] 6.4019321 4.3553425
[27,] -4.4811297 6.4019321
[28,] 6.1006296 -4.4811297
[29,] -3.3319566 6.1006296
[30,] 4.5207505 -3.3319566
[31,] -1.5444450 4.5207505
[32,] 2.7982492 -1.5444450
[33,] -8.9998609 2.7982492
[34,] 0.3877517 -8.9998609
[35,] 1.8130021 0.3877517
[36,] 9.2452093 1.8130021
[37,] 7.4768812 9.2452093
[38,] 0.5815520 7.4768812
[39,] 10.1441494 0.5815520
[40,] 7.7650979 10.1441494
[41,] 14.3369661 7.7650979
[42,] 13.2820224 14.3369661
[43,] 29.0089746 13.2820224
[44,] 10.7431419 29.0089746
[45,] 15.9890284 10.7431419
[46,] 1.8879596 15.9890284
[47,] 10.0647317 1.8879596
[48,] 5.3230073 10.0647317
[49,] 11.1011358 5.3230073
[50,] 0.1865522 11.1011358
[51,] 13.1371098 0.1865522
[52,] -14.6337651 13.1371098
[53,] -30.7868939 -14.6337651
[54,] -15.9176604 -30.7868939
[55,] -25.6838805 -15.9176604
[56,] -28.1259772 -25.6838805
[57,] -16.3980233 -28.1259772
[58,] -7.6813609 -16.3980233
[59,] -7.8683188 -7.6813609
[60,] -6.3156327 -7.8683188
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -19.0398810 -14.1708420
2 -13.1528185 -19.0398810
3 -14.6401740 -13.1528185
4 -10.8731540 -14.6401740
5 4.7160542 -10.8731540
6 -12.8267221 4.7160542
7 -9.8038500 -12.8267221
8 0.5101946 -9.8038500
9 1.8062297 0.5101946
10 -6.3822868 1.8062297
11 -7.3997630 -6.3822868
12 -1.0597898 -7.3997630
13 -3.8934786 -1.0597898
14 5.9827821 -3.8934786
15 -4.1599554 5.9827821
16 11.6411916 -4.1599554
17 15.0658302 11.6411916
18 10.9416096 15.0658302
19 8.0232010 10.9416096
20 14.0743916 8.0232010
21 7.6026260 14.0743916
22 11.7879363 7.6026260
23 3.3903480 11.7879363
24 6.9780478 3.3903480
25 4.3553425 6.9780478
26 6.4019321 4.3553425
27 -4.4811297 6.4019321
28 6.1006296 -4.4811297
29 -3.3319566 6.1006296
30 4.5207505 -3.3319566
31 -1.5444450 4.5207505
32 2.7982492 -1.5444450
33 -8.9998609 2.7982492
34 0.3877517 -8.9998609
35 1.8130021 0.3877517
36 9.2452093 1.8130021
37 7.4768812 9.2452093
38 0.5815520 7.4768812
39 10.1441494 0.5815520
40 7.7650979 10.1441494
41 14.3369661 7.7650979
42 13.2820224 14.3369661
43 29.0089746 13.2820224
44 10.7431419 29.0089746
45 15.9890284 10.7431419
46 1.8879596 15.9890284
47 10.0647317 1.8879596
48 5.3230073 10.0647317
49 11.1011358 5.3230073
50 0.1865522 11.1011358
51 13.1371098 0.1865522
52 -14.6337651 13.1371098
53 -30.7868939 -14.6337651
54 -15.9176604 -30.7868939
55 -25.6838805 -15.9176604
56 -28.1259772 -25.6838805
57 -16.3980233 -28.1259772
58 -7.6813609 -16.3980233
59 -7.8683188 -7.6813609
60 -6.3156327 -7.8683188
> 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/778061259087225.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/80dsn1259087225.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/9f31h1259087225.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/103q7m1259087225.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/11czmf1259087225.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/12381c1259087225.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/133mr61259087225.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/14p3xy1259087225.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/15x45l1259087225.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/1658ec1259087225.tab")
+ }
>
> system("convert tmp/111u91259087225.ps tmp/111u91259087225.png")
> system("convert tmp/2ro651259087225.ps tmp/2ro651259087225.png")
> system("convert tmp/3kksh1259087225.ps tmp/3kksh1259087225.png")
> system("convert tmp/4qn0e1259087225.ps tmp/4qn0e1259087225.png")
> system("convert tmp/5cy6x1259087225.ps tmp/5cy6x1259087225.png")
> system("convert tmp/6bffd1259087225.ps tmp/6bffd1259087225.png")
> system("convert tmp/778061259087225.ps tmp/778061259087225.png")
> system("convert tmp/80dsn1259087225.ps tmp/80dsn1259087225.png")
> system("convert tmp/9f31h1259087225.ps tmp/9f31h1259087225.png")
> system("convert tmp/103q7m1259087225.ps tmp/103q7m1259087225.png")
>
>
> proc.time()
user system elapsed
2.439 1.583 3.090