R version 2.8.0 (2008-10-20)
Copyright (C) 2008 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
Natural language support but running in an English locale
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> x <- array(list(37,1,30,1,47,0,35,1,30,1,43,0,82,0,40,0,47,0,19,1,52,0,136,0,80,0,42,0,54,0,66,0,81,0,63,0,137,0,72,0,107,0,58,0,36,1,52,0,79,0,77,0,54,0,84,0,48,0,96,0,83,0,66,0,61,0,53,0,30,1,74,0,69,0,59,0,42,0,65,0,70,0,100,0,63,0,105,0,82,0,81,0,75,0,102,0,121,0,98,0,76,0,77,0,63,0,37,1,35,1,23,1,40,0,29,1,37,1,51,0,20,1,28,1,13,1,22,1,25,1,13,1,16,1,13,1,16,1,17,1,9,1,17,1,25,1,14,1,8,1,7,1,10,1,7,1,10,1,3,1),dim=c(2,80),dimnames=list(c('SKIA','x'),1:80))
> y <- array(NA,dim=c(2,80),dimnames=list(c('SKIA','x'),1:80))
> for (i in 1:dim(x)[1])
+ {
+ for (j in 1:dim(x)[2])
+ {
+ y[i,j] <- as.numeric(x[i,j])
+ }
+ }
> par3 = '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
SKIA x M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 37 1 1 0 0 0 0 0 0 0 0 0 0 1
2 30 1 0 1 0 0 0 0 0 0 0 0 0 2
3 47 0 0 0 1 0 0 0 0 0 0 0 0 3
4 35 1 0 0 0 1 0 0 0 0 0 0 0 4
5 30 1 0 0 0 0 1 0 0 0 0 0 0 5
6 43 0 0 0 0 0 0 1 0 0 0 0 0 6
7 82 0 0 0 0 0 0 0 1 0 0 0 0 7
8 40 0 0 0 0 0 0 0 0 1 0 0 0 8
9 47 0 0 0 0 0 0 0 0 0 1 0 0 9
10 19 1 0 0 0 0 0 0 0 0 0 1 0 10
11 52 0 0 0 0 0 0 0 0 0 0 0 1 11
12 136 0 0 0 0 0 0 0 0 0 0 0 0 12
13 80 0 1 0 0 0 0 0 0 0 0 0 0 13
14 42 0 0 1 0 0 0 0 0 0 0 0 0 14
15 54 0 0 0 1 0 0 0 0 0 0 0 0 15
16 66 0 0 0 0 1 0 0 0 0 0 0 0 16
17 81 0 0 0 0 0 1 0 0 0 0 0 0 17
18 63 0 0 0 0 0 0 1 0 0 0 0 0 18
19 137 0 0 0 0 0 0 0 1 0 0 0 0 19
20 72 0 0 0 0 0 0 0 0 1 0 0 0 20
21 107 0 0 0 0 0 0 0 0 0 1 0 0 21
22 58 0 0 0 0 0 0 0 0 0 0 1 0 22
23 36 1 0 0 0 0 0 0 0 0 0 0 1 23
24 52 0 0 0 0 0 0 0 0 0 0 0 0 24
25 79 0 1 0 0 0 0 0 0 0 0 0 0 25
26 77 0 0 1 0 0 0 0 0 0 0 0 0 26
27 54 0 0 0 1 0 0 0 0 0 0 0 0 27
28 84 0 0 0 0 1 0 0 0 0 0 0 0 28
29 48 0 0 0 0 0 1 0 0 0 0 0 0 29
30 96 0 0 0 0 0 0 1 0 0 0 0 0 30
31 83 0 0 0 0 0 0 0 1 0 0 0 0 31
32 66 0 0 0 0 0 0 0 0 1 0 0 0 32
33 61 0 0 0 0 0 0 0 0 0 1 0 0 33
34 53 0 0 0 0 0 0 0 0 0 0 1 0 34
35 30 1 0 0 0 0 0 0 0 0 0 0 1 35
36 74 0 0 0 0 0 0 0 0 0 0 0 0 36
37 69 0 1 0 0 0 0 0 0 0 0 0 0 37
38 59 0 0 1 0 0 0 0 0 0 0 0 0 38
39 42 0 0 0 1 0 0 0 0 0 0 0 0 39
40 65 0 0 0 0 1 0 0 0 0 0 0 0 40
41 70 0 0 0 0 0 1 0 0 0 0 0 0 41
42 100 0 0 0 0 0 0 1 0 0 0 0 0 42
43 63 0 0 0 0 0 0 0 1 0 0 0 0 43
44 105 0 0 0 0 0 0 0 0 1 0 0 0 44
45 82 0 0 0 0 0 0 0 0 0 1 0 0 45
46 81 0 0 0 0 0 0 0 0 0 0 1 0 46
47 75 0 0 0 0 0 0 0 0 0 0 0 1 47
48 102 0 0 0 0 0 0 0 0 0 0 0 0 48
49 121 0 1 0 0 0 0 0 0 0 0 0 0 49
50 98 0 0 1 0 0 0 0 0 0 0 0 0 50
51 76 0 0 0 1 0 0 0 0 0 0 0 0 51
52 77 0 0 0 0 1 0 0 0 0 0 0 0 52
53 63 0 0 0 0 0 1 0 0 0 0 0 0 53
54 37 1 0 0 0 0 0 1 0 0 0 0 0 54
55 35 1 0 0 0 0 0 0 1 0 0 0 0 55
56 23 1 0 0 0 0 0 0 0 1 0 0 0 56
57 40 0 0 0 0 0 0 0 0 0 1 0 0 57
58 29 1 0 0 0 0 0 0 0 0 0 1 0 58
59 37 1 0 0 0 0 0 0 0 0 0 0 1 59
60 51 0 0 0 0 0 0 0 0 0 0 0 0 60
61 20 1 1 0 0 0 0 0 0 0 0 0 0 61
62 28 1 0 1 0 0 0 0 0 0 0 0 0 62
63 13 1 0 0 1 0 0 0 0 0 0 0 0 63
64 22 1 0 0 0 1 0 0 0 0 0 0 0 64
65 25 1 0 0 0 0 1 0 0 0 0 0 0 65
66 13 1 0 0 0 0 0 1 0 0 0 0 0 66
67 16 1 0 0 0 0 0 0 1 0 0 0 0 67
68 13 1 0 0 0 0 0 0 0 1 0 0 0 68
69 16 1 0 0 0 0 0 0 0 0 1 0 0 69
70 17 1 0 0 0 0 0 0 0 0 0 1 0 70
71 9 1 0 0 0 0 0 0 0 0 0 0 1 71
72 17 1 0 0 0 0 0 0 0 0 0 0 0 72
73 25 1 1 0 0 0 0 0 0 0 0 0 0 73
74 14 1 0 1 0 0 0 0 0 0 0 0 0 74
75 8 1 0 0 1 0 0 0 0 0 0 0 0 75
76 7 1 0 0 0 1 0 0 0 0 0 0 0 76
77 10 1 0 0 0 0 1 0 0 0 0 0 0 77
78 7 1 0 0 0 0 0 1 0 0 0 0 0 78
79 10 1 0 0 0 0 0 0 1 0 0 0 0 79
80 3 1 0 0 0 0 0 0 0 1 0 0 0 80
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) x M1 M2 M3 M4
84.3432 -49.4382 2.0310 -9.7284 -24.4076 -8.3902
M5 M6 M7 M8 M9 M10
-12.4353 -7.7662 1.9029 -12.8565 -13.4598 -12.8827
M11 t
-7.5453 -0.0977
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-32.991 -11.522 -2.382 8.877 52.829
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 84.3432 9.0787 9.290 1.33e-13 ***
x -49.4382 5.5031 -8.984 4.65e-13 ***
M1 2.0310 11.0270 0.184 0.8544
M2 -9.7284 11.0137 -0.883 0.3803
M3 -24.4076 10.9131 -2.237 0.0287 *
M4 -8.3902 10.9906 -0.763 0.4479
M5 -12.4353 10.9808 -1.132 0.2615
M6 -7.7662 10.9721 -0.708 0.4816
M7 1.9029 10.9646 0.174 0.8627
M8 -12.8565 10.9582 -1.173 0.2449
M9 -13.4598 11.2927 -1.192 0.2376
M10 -12.8827 11.4577 -1.124 0.2649
M11 -7.5453 11.6333 -0.649 0.5189
t -0.0977 0.1133 -0.863 0.3914
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 19.55 on 66 degrees of freedom
Multiple R-squared: 0.6866, Adjusted R-squared: 0.6248
F-statistic: 11.12 on 13 and 66 DF, p-value: 4.189e-12
> 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.2800873 5.601747e-01 7.199127e-01
[2,] 0.1758596 3.517193e-01 8.241404e-01
[3,] 0.4386295 8.772591e-01 5.613705e-01
[4,] 0.3125351 6.250703e-01 6.874649e-01
[5,] 0.4381612 8.763224e-01 5.618388e-01
[6,] 0.3413552 6.827104e-01 6.586448e-01
[7,] 0.4847382 9.694764e-01 5.152618e-01
[8,] 0.9940526 1.189483e-02 5.947416e-03
[9,] 0.9901421 1.971575e-02 9.857875e-03
[10,] 0.9837362 3.252755e-02 1.626377e-02
[11,] 0.9771851 4.562978e-02 2.281489e-02
[12,] 0.9638807 7.223855e-02 3.611928e-02
[13,] 0.9770166 4.596670e-02 2.298335e-02
[14,] 0.9744512 5.109764e-02 2.554882e-02
[15,] 0.9794374 4.112528e-02 2.056264e-02
[16,] 0.9716042 5.679152e-02 2.839576e-02
[17,] 0.9669954 6.600925e-02 3.300463e-02
[18,] 0.9700642 5.987151e-02 2.993575e-02
[19,] 0.9558237 8.835251e-02 4.417626e-02
[20,] 0.9504154 9.916914e-02 4.958457e-02
[21,] 0.9597892 8.042164e-02 4.021082e-02
[22,] 0.9735766 5.284672e-02 2.642336e-02
[23,] 0.9912824 1.743525e-02 8.717625e-03
[24,] 0.9927076 1.458474e-02 7.292371e-03
[25,] 0.9918757 1.624850e-02 8.124250e-03
[26,] 0.9910861 1.782778e-02 8.913892e-03
[27,] 0.9973055 5.389084e-03 2.694542e-03
[28,] 0.9986080 2.784088e-03 1.392044e-03
[29,] 0.9977467 4.506573e-03 2.253287e-03
[30,] 0.9963286 7.342758e-03 3.671379e-03
[31,] 0.9942906 1.141888e-02 5.709440e-03
[32,] 0.9960569 7.886204e-03 3.943102e-03
[33,] 0.9998657 2.686915e-04 1.343457e-04
[34,] 0.9999712 5.768834e-05 2.884417e-05
[35,] 0.9999823 3.546519e-05 1.773259e-05
[36,] 0.9999960 7.901277e-06 3.950638e-06
[37,] 0.9999921 1.586949e-05 7.934744e-06
[38,] 0.9999877 2.464988e-05 1.232494e-05
[39,] 0.9999717 5.667281e-05 2.833641e-05
[40,] 0.9998987 2.026725e-04 1.013362e-04
[41,] 0.9997963 4.073827e-04 2.036914e-04
[42,] 0.9992636 1.472871e-03 7.364354e-04
[43,] 0.9997654 4.691819e-04 2.345910e-04
[44,] 0.9991441 1.711820e-03 8.559100e-04
[45,] 0.9996620 6.760709e-04 3.380355e-04
[46,] 0.9983557 3.288628e-03 1.644314e-03
[47,] 0.9929004 1.419923e-02 7.099613e-03
> postscript(file="/var/www/html/freestat/rcomp/tmp/10vt91291143727.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> plot(x[,1], type='l', main='Actuals and Interpolation', ylab='value of Actuals and Interpolation (dots)', xlab='time or index')
> points(x[,1]-mysum$resid)
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/20vt91291143727.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> plot(mysum$resid, type='b', pch=19, main='Residuals', ylab='value of Residuals', xlab='time or index')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/30vt91291143727.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> hist(mysum$resid, main='Residual Histogram', xlab='values of Residuals')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/4a4au1291143727.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> densityplot(~mysum$resid,col='black',main='Residual Density Plot', xlab='values of Residuals')
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/5a4au1291143727.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> qqnorm(mysum$resid, main='Residual Normal Q-Q Plot')
> qqline(mysum$resid)
> grid()
> dev.off()
null device
1
> (myerror <- as.ts(mysum$resid))
Time Series:
Start = 1
End = 80
Frequency = 1
1 2 3 4 5 6
0.1617020 5.0188449 -12.6424577 8.8759877 8.0188449 -32.9907701
7 8 9 10 11 12
-3.5621986 -30.7050558 -23.0040824 -2.0452915 -23.7231756 52.8292509
13 14 15 16 17 18
-5.1040638 -31.2469209 -4.4700372 -8.3897781 10.7530791 -11.8183495
19 20 21 22 23 24
52.6102219 2.4673648 38.1683381 -11.3110573 10.8874313 -29.9983285
25 26 27 28 29 30
-4.9316433 4.9254996 -3.2976166 10.7826425 -21.0745004 22.3540710
31 32 33 34 35 36
-0.2173575 -2.3602147 -6.6592413 -15.1386368 6.0598518 -6.8259080
37 38 39 40 41 42
-13.7592227 -11.9020799 -14.1251961 -7.0449370 2.0979201 27.5264916
43 44 45 46 47 48
-19.0449370 37.8122058 15.5131792 14.0337838 2.7940860 22.3465125
49 50 51 52 53 54
39.4131978 28.2703407 21.0472244 6.1274835 -3.7296593 15.1370985
55 56 57 58 59 60
3.5656699 6.4228127 -25.3144002 12.6443907 15.4046929 -27.4810669
61 62 63 64 65 66
-10.9761953 8.8809476 8.6578313 1.7380904 8.8809476 -7.6904810
67 68 69 70 71 72
-14.2619096 -2.4047667 1.2962066 1.8168112 -11.4228865 -10.8704600
73 74 75 76 77 78
-4.8037747 -3.9466319 4.8302519 -12.0894890 -4.9466319 -12.5180605
79 80
-19.0894890 -11.2323462
> postscript(file="/var/www/html/freestat/rcomp/tmp/6a4au1291143727.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> dum <- cbind(lag(myerror,k=1),myerror)
> dum
Time Series:
Start = 0
End = 80
Frequency = 1
lag(myerror, k = 1) myerror
0 0.1617020 NA
1 5.0188449 0.1617020
2 -12.6424577 5.0188449
3 8.8759877 -12.6424577
4 8.0188449 8.8759877
5 -32.9907701 8.0188449
6 -3.5621986 -32.9907701
7 -30.7050558 -3.5621986
8 -23.0040824 -30.7050558
9 -2.0452915 -23.0040824
10 -23.7231756 -2.0452915
11 52.8292509 -23.7231756
12 -5.1040638 52.8292509
13 -31.2469209 -5.1040638
14 -4.4700372 -31.2469209
15 -8.3897781 -4.4700372
16 10.7530791 -8.3897781
17 -11.8183495 10.7530791
18 52.6102219 -11.8183495
19 2.4673648 52.6102219
20 38.1683381 2.4673648
21 -11.3110573 38.1683381
22 10.8874313 -11.3110573
23 -29.9983285 10.8874313
24 -4.9316433 -29.9983285
25 4.9254996 -4.9316433
26 -3.2976166 4.9254996
27 10.7826425 -3.2976166
28 -21.0745004 10.7826425
29 22.3540710 -21.0745004
30 -0.2173575 22.3540710
31 -2.3602147 -0.2173575
32 -6.6592413 -2.3602147
33 -15.1386368 -6.6592413
34 6.0598518 -15.1386368
35 -6.8259080 6.0598518
36 -13.7592227 -6.8259080
37 -11.9020799 -13.7592227
38 -14.1251961 -11.9020799
39 -7.0449370 -14.1251961
40 2.0979201 -7.0449370
41 27.5264916 2.0979201
42 -19.0449370 27.5264916
43 37.8122058 -19.0449370
44 15.5131792 37.8122058
45 14.0337838 15.5131792
46 2.7940860 14.0337838
47 22.3465125 2.7940860
48 39.4131978 22.3465125
49 28.2703407 39.4131978
50 21.0472244 28.2703407
51 6.1274835 21.0472244
52 -3.7296593 6.1274835
53 15.1370985 -3.7296593
54 3.5656699 15.1370985
55 6.4228127 3.5656699
56 -25.3144002 6.4228127
57 12.6443907 -25.3144002
58 15.4046929 12.6443907
59 -27.4810669 15.4046929
60 -10.9761953 -27.4810669
61 8.8809476 -10.9761953
62 8.6578313 8.8809476
63 1.7380904 8.6578313
64 8.8809476 1.7380904
65 -7.6904810 8.8809476
66 -14.2619096 -7.6904810
67 -2.4047667 -14.2619096
68 1.2962066 -2.4047667
69 1.8168112 1.2962066
70 -11.4228865 1.8168112
71 -10.8704600 -11.4228865
72 -4.8037747 -10.8704600
73 -3.9466319 -4.8037747
74 4.8302519 -3.9466319
75 -12.0894890 4.8302519
76 -4.9466319 -12.0894890
77 -12.5180605 -4.9466319
78 -19.0894890 -12.5180605
79 -11.2323462 -19.0894890
80 NA -11.2323462
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 5.0188449 0.1617020
[2,] -12.6424577 5.0188449
[3,] 8.8759877 -12.6424577
[4,] 8.0188449 8.8759877
[5,] -32.9907701 8.0188449
[6,] -3.5621986 -32.9907701
[7,] -30.7050558 -3.5621986
[8,] -23.0040824 -30.7050558
[9,] -2.0452915 -23.0040824
[10,] -23.7231756 -2.0452915
[11,] 52.8292509 -23.7231756
[12,] -5.1040638 52.8292509
[13,] -31.2469209 -5.1040638
[14,] -4.4700372 -31.2469209
[15,] -8.3897781 -4.4700372
[16,] 10.7530791 -8.3897781
[17,] -11.8183495 10.7530791
[18,] 52.6102219 -11.8183495
[19,] 2.4673648 52.6102219
[20,] 38.1683381 2.4673648
[21,] -11.3110573 38.1683381
[22,] 10.8874313 -11.3110573
[23,] -29.9983285 10.8874313
[24,] -4.9316433 -29.9983285
[25,] 4.9254996 -4.9316433
[26,] -3.2976166 4.9254996
[27,] 10.7826425 -3.2976166
[28,] -21.0745004 10.7826425
[29,] 22.3540710 -21.0745004
[30,] -0.2173575 22.3540710
[31,] -2.3602147 -0.2173575
[32,] -6.6592413 -2.3602147
[33,] -15.1386368 -6.6592413
[34,] 6.0598518 -15.1386368
[35,] -6.8259080 6.0598518
[36,] -13.7592227 -6.8259080
[37,] -11.9020799 -13.7592227
[38,] -14.1251961 -11.9020799
[39,] -7.0449370 -14.1251961
[40,] 2.0979201 -7.0449370
[41,] 27.5264916 2.0979201
[42,] -19.0449370 27.5264916
[43,] 37.8122058 -19.0449370
[44,] 15.5131792 37.8122058
[45,] 14.0337838 15.5131792
[46,] 2.7940860 14.0337838
[47,] 22.3465125 2.7940860
[48,] 39.4131978 22.3465125
[49,] 28.2703407 39.4131978
[50,] 21.0472244 28.2703407
[51,] 6.1274835 21.0472244
[52,] -3.7296593 6.1274835
[53,] 15.1370985 -3.7296593
[54,] 3.5656699 15.1370985
[55,] 6.4228127 3.5656699
[56,] -25.3144002 6.4228127
[57,] 12.6443907 -25.3144002
[58,] 15.4046929 12.6443907
[59,] -27.4810669 15.4046929
[60,] -10.9761953 -27.4810669
[61,] 8.8809476 -10.9761953
[62,] 8.6578313 8.8809476
[63,] 1.7380904 8.6578313
[64,] 8.8809476 1.7380904
[65,] -7.6904810 8.8809476
[66,] -14.2619096 -7.6904810
[67,] -2.4047667 -14.2619096
[68,] 1.2962066 -2.4047667
[69,] 1.8168112 1.2962066
[70,] -11.4228865 1.8168112
[71,] -10.8704600 -11.4228865
[72,] -4.8037747 -10.8704600
[73,] -3.9466319 -4.8037747
[74,] 4.8302519 -3.9466319
[75,] -12.0894890 4.8302519
[76,] -4.9466319 -12.0894890
[77,] -12.5180605 -4.9466319
[78,] -19.0894890 -12.5180605
[79,] -11.2323462 -19.0894890
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 5.0188449 0.1617020
2 -12.6424577 5.0188449
3 8.8759877 -12.6424577
4 8.0188449 8.8759877
5 -32.9907701 8.0188449
6 -3.5621986 -32.9907701
7 -30.7050558 -3.5621986
8 -23.0040824 -30.7050558
9 -2.0452915 -23.0040824
10 -23.7231756 -2.0452915
11 52.8292509 -23.7231756
12 -5.1040638 52.8292509
13 -31.2469209 -5.1040638
14 -4.4700372 -31.2469209
15 -8.3897781 -4.4700372
16 10.7530791 -8.3897781
17 -11.8183495 10.7530791
18 52.6102219 -11.8183495
19 2.4673648 52.6102219
20 38.1683381 2.4673648
21 -11.3110573 38.1683381
22 10.8874313 -11.3110573
23 -29.9983285 10.8874313
24 -4.9316433 -29.9983285
25 4.9254996 -4.9316433
26 -3.2976166 4.9254996
27 10.7826425 -3.2976166
28 -21.0745004 10.7826425
29 22.3540710 -21.0745004
30 -0.2173575 22.3540710
31 -2.3602147 -0.2173575
32 -6.6592413 -2.3602147
33 -15.1386368 -6.6592413
34 6.0598518 -15.1386368
35 -6.8259080 6.0598518
36 -13.7592227 -6.8259080
37 -11.9020799 -13.7592227
38 -14.1251961 -11.9020799
39 -7.0449370 -14.1251961
40 2.0979201 -7.0449370
41 27.5264916 2.0979201
42 -19.0449370 27.5264916
43 37.8122058 -19.0449370
44 15.5131792 37.8122058
45 14.0337838 15.5131792
46 2.7940860 14.0337838
47 22.3465125 2.7940860
48 39.4131978 22.3465125
49 28.2703407 39.4131978
50 21.0472244 28.2703407
51 6.1274835 21.0472244
52 -3.7296593 6.1274835
53 15.1370985 -3.7296593
54 3.5656699 15.1370985
55 6.4228127 3.5656699
56 -25.3144002 6.4228127
57 12.6443907 -25.3144002
58 15.4046929 12.6443907
59 -27.4810669 15.4046929
60 -10.9761953 -27.4810669
61 8.8809476 -10.9761953
62 8.6578313 8.8809476
63 1.7380904 8.6578313
64 8.8809476 1.7380904
65 -7.6904810 8.8809476
66 -14.2619096 -7.6904810
67 -2.4047667 -14.2619096
68 1.2962066 -2.4047667
69 1.8168112 1.2962066
70 -11.4228865 1.8168112
71 -10.8704600 -11.4228865
72 -4.8037747 -10.8704600
73 -3.9466319 -4.8037747
74 4.8302519 -3.9466319
75 -12.0894890 4.8302519
76 -4.9466319 -12.0894890
77 -12.5180605 -4.9466319
78 -19.0894890 -12.5180605
79 -11.2323462 -19.0894890
> 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/7leax1291143727.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> acf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Autocorrelation Function')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/8wn901291143727.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> pacf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Partial Autocorrelation Function')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/9wn901291143727.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> opar <- par(mfrow = c(2,2), oma = c(0, 0, 1.1, 0))
> plot(mylm, las = 1, sub='Residual Diagnostics')
> par(opar)
> dev.off()
null device
1
> if (n > n25) {
+ postscript(file="/var/www/html/freestat/rcomp/tmp/10pw831291143727.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
+ plot(kp3:nmkm3,gqarr[,2], main='Goldfeld-Quandt test',ylab='2-sided p-value',xlab='breakpoint')
+ grid()
+ dev.off()
+ }
null device
1
>
> #Note: the /var/www/html/freestat/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/html/freestat/rcomp/createtable")
>
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Estimated Regression Equation', 1, TRUE)
> a<-table.row.end(a)
> myeq <- colnames(x)[1]
> myeq <- paste(myeq, '[t] = ', sep='')
> for (i in 1:k){
+ if (mysum$coefficients[i,1] > 0) myeq <- paste(myeq, '+', '')
+ myeq <- paste(myeq, mysum$coefficients[i,1], sep=' ')
+ if (rownames(mysum$coefficients)[i] != '(Intercept)') {
+ myeq <- paste(myeq, rownames(mysum$coefficients)[i], sep='')
+ if (rownames(mysum$coefficients)[i] != 't') myeq <- paste(myeq, '[t]', sep='')
+ }
+ }
> myeq <- paste(myeq, ' + e[t]')
> a<-table.row.start(a)
> a<-table.element(a, myeq)
> a<-table.row.end(a)
> a<-table.end(a)
> table.save(a,file="/var/www/html/freestat/rcomp/tmp/11sfpr1291143727.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/12df5e1291143727.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/13r7l51291143727.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/145i4o1291143728.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/15r02c1291143728.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/16uj1i1291143728.tab")
+ }
>
> try(system("convert tmp/10vt91291143727.ps tmp/10vt91291143727.png",intern=TRUE))
character(0)
> try(system("convert tmp/20vt91291143727.ps tmp/20vt91291143727.png",intern=TRUE))
character(0)
> try(system("convert tmp/30vt91291143727.ps tmp/30vt91291143727.png",intern=TRUE))
character(0)
> try(system("convert tmp/4a4au1291143727.ps tmp/4a4au1291143727.png",intern=TRUE))
character(0)
> try(system("convert tmp/5a4au1291143727.ps tmp/5a4au1291143727.png",intern=TRUE))
character(0)
> try(system("convert tmp/6a4au1291143727.ps tmp/6a4au1291143727.png",intern=TRUE))
character(0)
> try(system("convert tmp/7leax1291143727.ps tmp/7leax1291143727.png",intern=TRUE))
character(0)
> try(system("convert tmp/8wn901291143727.ps tmp/8wn901291143727.png",intern=TRUE))
character(0)
> try(system("convert tmp/9wn901291143727.ps tmp/9wn901291143727.png",intern=TRUE))
character(0)
> try(system("convert tmp/10pw831291143727.ps tmp/10pw831291143727.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
4.026 2.456 4.483