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(16571771.60
+ ,17972385.83
+ ,17535166.38
+ ,16198106.13
+ ,17487530.67
+ ,13768040.14
+ ,16198892.67
+ ,16896235.55
+ ,16571771.60
+ ,17535166.38
+ ,16198106.13
+ ,17487530.67
+ ,16554237.93
+ ,16697955.94
+ ,16198892.67
+ ,16571771.60
+ ,17535166.38
+ ,16198106.13
+ ,19554176.37
+ ,19691579.52
+ ,16554237.93
+ ,16198892.67
+ ,16571771.60
+ ,17535166.38
+ ,15903762.33
+ ,15930700.75
+ ,19554176.37
+ ,16554237.93
+ ,16198892.67
+ ,16571771.60
+ ,18003781.65
+ ,17444615.98
+ ,15903762.33
+ ,19554176.37
+ ,16554237.93
+ ,16198892.67
+ ,18329610.38
+ ,17699369.88
+ ,18003781.65
+ ,15903762.33
+ ,19554176.37
+ ,16554237.93
+ ,16260733.42
+ ,15189796.81
+ ,18329610.38
+ ,18003781.65
+ ,15903762.33
+ ,19554176.37
+ ,14851949.20
+ ,15672722.75
+ ,16260733.42
+ ,18329610.38
+ ,18003781.65
+ ,15903762.33
+ ,18174068.44
+ ,17180794.3
+ ,14851949.20
+ ,16260733.42
+ ,18329610.38
+ ,18003781.65
+ ,18406552.23
+ ,17664893.45
+ ,18174068.44
+ ,14851949.20
+ ,16260733.42
+ ,18329610.38
+ ,18466459.42
+ ,17862884.98
+ ,18406552.23
+ ,18174068.44
+ ,14851949.20
+ ,16260733.42
+ ,16016524.60
+ ,16162288.88
+ ,18466459.42
+ ,18406552.23
+ ,18174068.44
+ ,14851949.20
+ ,17428458.32
+ ,17463628.82
+ ,16016524.60
+ ,18466459.42
+ ,18406552.23
+ ,18174068.44
+ ,17167191.42
+ ,16772112.17
+ ,17428458.32
+ ,16016524.60
+ ,18466459.42
+ ,18406552.23
+ ,19629987.60
+ ,19106861.48
+ ,17167191.42
+ ,17428458.32
+ ,16016524.60
+ ,18466459.42
+ ,17183629.01
+ ,16721314.25
+ ,19629987.60
+ ,17167191.42
+ ,17428458.32
+ ,16016524.60
+ ,18344657.85
+ ,18161267.85
+ ,17183629.01
+ ,19629987.60
+ ,17167191.42
+ ,17428458.32
+ ,19301440.71
+ ,18509941.2
+ ,18344657.85
+ ,17183629.01
+ ,19629987.60
+ ,17167191.42
+ ,18147463.68
+ ,17802737.97
+ ,19301440.71
+ ,18344657.85
+ ,17183629.01
+ ,19629987.60
+ ,16192909.22
+ ,16409869.75
+ ,18147463.68
+ ,19301440.71
+ ,18344657.85
+ ,17183629.01
+ ,18374420.60
+ ,17967742.04
+ ,16192909.22
+ ,18147463.68
+ ,19301440.71
+ ,18344657.85
+ ,20515191.95
+ ,20286602.27
+ ,18374420.60
+ ,16192909.22
+ ,18147463.68
+ ,19301440.71
+ ,18957217.20
+ ,19537280.81
+ ,20515191.95
+ ,18374420.60
+ ,16192909.22
+ ,18147463.68
+ ,16471529.53
+ ,18021889.62
+ ,18957217.20
+ ,20515191.95
+ ,18374420.60
+ ,16192909.22
+ ,18746813.27
+ ,20194317.23
+ ,16471529.53
+ ,18957217.20
+ ,20515191.95
+ ,18374420.60
+ ,19009453.59
+ ,19049596.62
+ ,18746813.27
+ ,16471529.53
+ ,18957217.20
+ ,20515191.95
+ ,19211178.55
+ ,20244720.94
+ ,19009453.59
+ ,18746813.27
+ ,16471529.53
+ ,18957217.20
+ ,20547653.75
+ ,21473302.24
+ ,19211178.55
+ ,19009453.59
+ ,18746813.27
+ ,16471529.53
+ ,19325754.03
+ ,19673603.19
+ ,20547653.75
+ ,19211178.55
+ ,19009453.59
+ ,18746813.27
+ ,20605542.58
+ ,21053177.29
+ ,19325754.03
+ ,20547653.75
+ ,19211178.55
+ ,19009453.59
+ ,20056915.06
+ ,20159479.84
+ ,20605542.58
+ ,19325754.03
+ ,20547653.75
+ ,19211178.55
+ ,16141449.72
+ ,18203628.31
+ ,20056915.06
+ ,20605542.58
+ ,19325754.03
+ ,20547653.75
+ ,20359793.22
+ ,21289464.94
+ ,16141449.72
+ ,20056915.06
+ ,20605542.58
+ ,19325754.03
+ ,19711553.27
+ ,20432335.71
+ ,20359793.22
+ ,16141449.72
+ ,20056915.06
+ ,20605542.58
+ ,15638580.70
+ ,17180395.07
+ ,19711553.27
+ ,20359793.22
+ ,16141449.72
+ ,20056915.06
+ ,14384486.00
+ ,15816786.32
+ ,15638580.70
+ ,19711553.27
+ ,20359793.22
+ ,16141449.72
+ ,13855616.12
+ ,15071819.75
+ ,14384486.00
+ ,15638580.70
+ ,19711553.27
+ ,20359793.22
+ ,14308336.46
+ ,14521120.61
+ ,13855616.12
+ ,14384486.00
+ ,15638580.70
+ ,19711553.27
+ ,15290621.44
+ ,15668789.39
+ ,14308336.46
+ ,13855616.12
+ ,14384486.00
+ ,15638580.70
+ ,14423755.53
+ ,14346884.11
+ ,15290621.44
+ ,14308336.46
+ ,13855616.12
+ ,14384486.00
+ ,13779681.49
+ ,13881008.13
+ ,14423755.53
+ ,15290621.44
+ ,14308336.46
+ ,13855616.12
+ ,15686348.94
+ ,15465943.69
+ ,13779681.49
+ ,14423755.53
+ ,15290621.44
+ ,14308336.46
+ ,14733828.17
+ ,14238232.92
+ ,15686348.94
+ ,13779681.49
+ ,14423755.53
+ ,15290621.44
+ ,12522497.94
+ ,13557713.21
+ ,14733828.17
+ ,15686348.94
+ ,13779681.49
+ ,14423755.53
+ ,16189383.57
+ ,16127590.29
+ ,12522497.94
+ ,14733828.17
+ ,15686348.94
+ ,13779681.49
+ ,16059123.25
+ ,16793894.2
+ ,16189383.57
+ ,12522497.94
+ ,14733828.17
+ ,15686348.94
+ ,16007123.26
+ ,16014007.43
+ ,16059123.25
+ ,16189383.57
+ ,12522497.94
+ ,14733828.17
+ ,15806842.33
+ ,16867867.15
+ ,16007123.26
+ ,16059123.25
+ ,16189383.57
+ ,12522497.94
+ ,15159951.13
+ ,16014583.21
+ ,15806842.33
+ ,16007123.26
+ ,16059123.25
+ ,16189383.57
+ ,15692144.17
+ ,15878594.85
+ ,15159951.13
+ ,15806842.33
+ ,16007123.26
+ ,16059123.25
+ ,18908869.11
+ ,18664899.14
+ ,15692144.17
+ ,15159951.13
+ ,15806842.33
+ ,16007123.26
+ ,16969881.42
+ ,17962530.06
+ ,18908869.11
+ ,15692144.17
+ ,15159951.13
+ ,15806842.33
+ ,16997477.78
+ ,17332692.2
+ ,16969881.42
+ ,18908869.11
+ ,15692144.17
+ ,15159951.13
+ ,19858875.65
+ ,19542066.35
+ ,16997477.78
+ ,16969881.42
+ ,18908869.11
+ ,15692144.17
+ ,17681170.13
+ ,17203555.19
+ ,19858875.65
+ ,16997477.78
+ ,16969881.42
+ ,18908869.11)
+ ,dim=c(6
+ ,56)
+ ,dimnames=list(c('Y'
+ ,'X'
+ ,'Y1'
+ ,'Y2'
+ ,'Y3'
+ ,'Y4')
+ ,1:56))
> y <- array(NA,dim=c(6,56),dimnames=list(c('Y','X','Y1','Y2','Y3','Y4'),1:56))
> 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
Y X Y1 Y2 Y3 Y4 M1 M2 M3 M4 M5 M6 M7
1 16571772 17972386 17535166 16198106 17487531 13768040 1 0 0 0 0 0 0
2 16198893 16896236 16571772 17535166 16198106 17487531 0 1 0 0 0 0 0
3 16554238 16697956 16198893 16571772 17535166 16198106 0 0 1 0 0 0 0
4 19554176 19691580 16554238 16198893 16571772 17535166 0 0 0 1 0 0 0
5 15903762 15930701 19554176 16554238 16198893 16571772 0 0 0 0 1 0 0
6 18003782 17444616 15903762 19554176 16554238 16198893 0 0 0 0 0 1 0
7 18329610 17699370 18003782 15903762 19554176 16554238 0 0 0 0 0 0 1
8 16260733 15189797 18329610 18003782 15903762 19554176 0 0 0 0 0 0 0
9 14851949 15672723 16260733 18329610 18003782 15903762 0 0 0 0 0 0 0
10 18174068 17180794 14851949 16260733 18329610 18003782 0 0 0 0 0 0 0
11 18406552 17664893 18174068 14851949 16260733 18329610 0 0 0 0 0 0 0
12 18466459 17862885 18406552 18174068 14851949 16260733 0 0 0 0 0 0 0
13 16016525 16162289 18466459 18406552 18174068 14851949 1 0 0 0 0 0 0
14 17428458 17463629 16016525 18466459 18406552 18174068 0 1 0 0 0 0 0
15 17167191 16772112 17428458 16016525 18466459 18406552 0 0 1 0 0 0 0
16 19629988 19106861 17167191 17428458 16016525 18466459 0 0 0 1 0 0 0
17 17183629 16721314 19629988 17167191 17428458 16016525 0 0 0 0 1 0 0
18 18344658 18161268 17183629 19629988 17167191 17428458 0 0 0 0 0 1 0
19 19301441 18509941 18344658 17183629 19629988 17167191 0 0 0 0 0 0 1
20 18147464 17802738 19301441 18344658 17183629 19629988 0 0 0 0 0 0 0
21 16192909 16409870 18147464 19301441 18344658 17183629 0 0 0 0 0 0 0
22 18374421 17967742 16192909 18147464 19301441 18344658 0 0 0 0 0 0 0
23 20515192 20286602 18374421 16192909 18147464 19301441 0 0 0 0 0 0 0
24 18957217 19537281 20515192 18374421 16192909 18147464 0 0 0 0 0 0 0
25 16471530 18021890 18957217 20515192 18374421 16192909 1 0 0 0 0 0 0
26 18746813 20194317 16471530 18957217 20515192 18374421 0 1 0 0 0 0 0
27 19009454 19049597 18746813 16471530 18957217 20515192 0 0 1 0 0 0 0
28 19211179 20244721 19009454 18746813 16471530 18957217 0 0 0 1 0 0 0
29 20547654 21473302 19211179 19009454 18746813 16471530 0 0 0 0 1 0 0
30 19325754 19673603 20547654 19211179 19009454 18746813 0 0 0 0 0 1 0
31 20605543 21053177 19325754 20547654 19211179 19009454 0 0 0 0 0 0 1
32 20056915 20159480 20605543 19325754 20547654 19211179 0 0 0 0 0 0 0
33 16141450 18203628 20056915 20605543 19325754 20547654 0 0 0 0 0 0 0
34 20359793 21289465 16141450 20056915 20605543 19325754 0 0 0 0 0 0 0
35 19711553 20432336 20359793 16141450 20056915 20605543 0 0 0 0 0 0 0
36 15638581 17180395 19711553 20359793 16141450 20056915 0 0 0 0 0 0 0
37 14384486 15816786 15638581 19711553 20359793 16141450 1 0 0 0 0 0 0
38 13855616 15071820 14384486 15638581 19711553 20359793 0 1 0 0 0 0 0
39 14308336 14521121 13855616 14384486 15638581 19711553 0 0 1 0 0 0 0
40 15290621 15668789 14308336 13855616 14384486 15638581 0 0 0 1 0 0 0
41 14423756 14346884 15290621 14308336 13855616 14384486 0 0 0 0 1 0 0
42 13779681 13881008 14423756 15290621 14308336 13855616 0 0 0 0 0 1 0
43 15686349 15465944 13779681 14423756 15290621 14308336 0 0 0 0 0 0 1
44 14733828 14238233 15686349 13779681 14423756 15290621 0 0 0 0 0 0 0
45 12522498 13557713 14733828 15686349 13779681 14423756 0 0 0 0 0 0 0
46 16189384 16127590 12522498 14733828 15686349 13779681 0 0 0 0 0 0 0
47 16059123 16793894 16189384 12522498 14733828 15686349 0 0 0 0 0 0 0
48 16007123 16014007 16059123 16189384 12522498 14733828 0 0 0 0 0 0 0
49 15806842 16867867 16007123 16059123 16189384 12522498 1 0 0 0 0 0 0
50 15159951 16014583 15806842 16007123 16059123 16189384 0 1 0 0 0 0 0
51 15692144 15878595 15159951 15806842 16007123 16059123 0 0 1 0 0 0 0
52 18908869 18664899 15692144 15159951 15806842 16007123 0 0 0 1 0 0 0
53 16969881 17962530 18908869 15692144 15159951 15806842 0 0 0 0 1 0 0
54 16997478 17332692 16969881 18908869 15692144 15159951 0 0 0 0 0 1 0
55 19858876 19542066 16997478 16969881 18908869 15692144 0 0 0 0 0 0 1
56 17681170 17203555 19858876 16997478 16969881 18908869 0 0 0 0 0 0 0
M8 M9 M10 M11 t
1 0 0 0 0 1
2 0 0 0 0 2
3 0 0 0 0 3
4 0 0 0 0 4
5 0 0 0 0 5
6 0 0 0 0 6
7 0 0 0 0 7
8 1 0 0 0 8
9 0 1 0 0 9
10 0 0 1 0 10
11 0 0 0 1 11
12 0 0 0 0 12
13 0 0 0 0 13
14 0 0 0 0 14
15 0 0 0 0 15
16 0 0 0 0 16
17 0 0 0 0 17
18 0 0 0 0 18
19 0 0 0 0 19
20 1 0 0 0 20
21 0 1 0 0 21
22 0 0 1 0 22
23 0 0 0 1 23
24 0 0 0 0 24
25 0 0 0 0 25
26 0 0 0 0 26
27 0 0 0 0 27
28 0 0 0 0 28
29 0 0 0 0 29
30 0 0 0 0 30
31 0 0 0 0 31
32 1 0 0 0 32
33 0 1 0 0 33
34 0 0 1 0 34
35 0 0 0 1 35
36 0 0 0 0 36
37 0 0 0 0 37
38 0 0 0 0 38
39 0 0 0 0 39
40 0 0 0 0 40
41 0 0 0 0 41
42 0 0 0 0 42
43 0 0 0 0 43
44 1 0 0 0 44
45 0 1 0 0 45
46 0 0 1 0 46
47 0 0 0 1 47
48 0 0 0 0 48
49 0 0 0 0 49
50 0 0 0 0 50
51 0 0 0 0 51
52 0 0 0 0 52
53 0 0 0 0 53
54 0 0 0 0 54
55 0 0 0 0 55
56 1 0 0 0 56
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X Y1 Y2 Y3 Y4
3.151e+06 8.625e-01 5.707e-02 -9.101e-03 3.989e-02 -1.198e-01
M1 M2 M3 M4 M5 M6
-1.281e+06 -4.903e+05 2.687e+05 3.998e+05 -2.017e+05 2.352e+05
M7 M8 M9 M10 M11 t
6.458e+05 8.228e+05 -9.885e+05 6.176e+05 4.390e+05 -1.761e+04
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-877749 -259402 -64440 346351 713058
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.151e+06 9.356e+05 3.368 0.001744 **
X 8.625e-01 7.859e-02 10.975 2.42e-13 ***
Y1 5.707e-02 8.025e-02 0.711 0.481321
Y2 -9.101e-03 7.935e-02 -0.115 0.909298
Y3 3.989e-02 7.823e-02 0.510 0.613008
Y4 -1.198e-01 7.454e-02 -1.607 0.116365
M1 -1.281e+06 5.091e+05 -2.517 0.016180 *
M2 -4.903e+05 4.808e+05 -1.020 0.314249
M3 2.687e+05 4.662e+05 0.577 0.567672
M4 3.998e+05 4.402e+05 0.908 0.369485
M5 -2.017e+05 3.861e+05 -0.522 0.604419
M6 2.352e+05 3.877e+05 0.607 0.547683
M7 6.458e+05 4.660e+05 1.386 0.173846
M8 8.228e+05 3.800e+05 2.165 0.036709 *
M9 -9.885e+05 4.435e+05 -2.229 0.031815 *
M10 6.176e+05 5.483e+05 1.126 0.267019
M11 4.390e+05 5.180e+05 0.848 0.402010
t -1.761e+04 4.518e+03 -3.899 0.000381 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 486000 on 38 degrees of freedom
Multiple R-squared: 0.9593, Adjusted R-squared: 0.9411
F-statistic: 52.71 on 17 and 38 DF, p-value: < 2.2e-16
> 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.3305387 0.6610773 0.6694613
[2,] 0.4738821 0.9477641 0.5261179
[3,] 0.5557129 0.8885742 0.4442871
[4,] 0.5420428 0.9159144 0.4579572
[5,] 0.5863171 0.8273658 0.4136829
[6,] 0.5566791 0.8866418 0.4433209
[7,] 0.5433760 0.9132481 0.4566240
[8,] 0.6041701 0.7916599 0.3958299
[9,] 0.5414989 0.9170023 0.4585011
[10,] 0.5206830 0.9586340 0.4793170
[11,] 0.4385472 0.8770944 0.5614528
[12,] 0.3585320 0.7170639 0.6414680
[13,] 0.2599881 0.5199762 0.7400119
[14,] 0.2059480 0.4118960 0.7940520
[15,] 0.4339922 0.8679845 0.5660078
> postscript(file="/var/www/html/rcomp/tmp/1eouq1290778168.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/2eouq1290778168.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/3pxta1290778168.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/4pxta1290778168.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/5pxta1290778168.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 = 56
Frequency = 1
1 2 3 4 5 6
-683405.991 -337410.474 -747747.286 -268247.707 -324413.885 227418.914
7 8 9 10 11 12
-289635.943 151957.740 -244353.058 488656.727 418988.201 590164.425
13 14 15 16 17 18
603223.234 648438.955 164727.429 633117.096 370634.085 212022.909
19 20 21 22 23 24
257045.108 -97859.522 713058.296 164350.491 519750.260 -97870.335
25 26 27 28 29 30
-190459.271 -258595.195 415969.927 -608210.939 -109831.837 -11252.797
31 32 33 34 35 36
-208968.725 -259515.312 -407440.176 -417992.914 -231982.145 -877749.447
37 38 39 40 41 42
-67614.913 -161722.783 128145.511 -461260.720 250036.215 -434512.183
43 44 45 46 47 48
-243868.920 -259364.892 -61265.061 -235014.305 -706756.316 385455.358
49 50 51 52 53 54
338256.940 109289.497 38904.419 704602.269 -186424.578 6323.156
55 56
485428.480 464781.987
> postscript(file="/var/www/html/rcomp/tmp/6iotd1290778168.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 = 56
Frequency = 1
lag(myerror, k = 1) myerror
0 -683405.991 NA
1 -337410.474 -683405.991
2 -747747.286 -337410.474
3 -268247.707 -747747.286
4 -324413.885 -268247.707
5 227418.914 -324413.885
6 -289635.943 227418.914
7 151957.740 -289635.943
8 -244353.058 151957.740
9 488656.727 -244353.058
10 418988.201 488656.727
11 590164.425 418988.201
12 603223.234 590164.425
13 648438.955 603223.234
14 164727.429 648438.955
15 633117.096 164727.429
16 370634.085 633117.096
17 212022.909 370634.085
18 257045.108 212022.909
19 -97859.522 257045.108
20 713058.296 -97859.522
21 164350.491 713058.296
22 519750.260 164350.491
23 -97870.335 519750.260
24 -190459.271 -97870.335
25 -258595.195 -190459.271
26 415969.927 -258595.195
27 -608210.939 415969.927
28 -109831.837 -608210.939
29 -11252.797 -109831.837
30 -208968.725 -11252.797
31 -259515.312 -208968.725
32 -407440.176 -259515.312
33 -417992.914 -407440.176
34 -231982.145 -417992.914
35 -877749.447 -231982.145
36 -67614.913 -877749.447
37 -161722.783 -67614.913
38 128145.511 -161722.783
39 -461260.720 128145.511
40 250036.215 -461260.720
41 -434512.183 250036.215
42 -243868.920 -434512.183
43 -259364.892 -243868.920
44 -61265.061 -259364.892
45 -235014.305 -61265.061
46 -706756.316 -235014.305
47 385455.358 -706756.316
48 338256.940 385455.358
49 109289.497 338256.940
50 38904.419 109289.497
51 704602.269 38904.419
52 -186424.578 704602.269
53 6323.156 -186424.578
54 485428.480 6323.156
55 464781.987 485428.480
56 NA 464781.987
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -337410.474 -683405.991
[2,] -747747.286 -337410.474
[3,] -268247.707 -747747.286
[4,] -324413.885 -268247.707
[5,] 227418.914 -324413.885
[6,] -289635.943 227418.914
[7,] 151957.740 -289635.943
[8,] -244353.058 151957.740
[9,] 488656.727 -244353.058
[10,] 418988.201 488656.727
[11,] 590164.425 418988.201
[12,] 603223.234 590164.425
[13,] 648438.955 603223.234
[14,] 164727.429 648438.955
[15,] 633117.096 164727.429
[16,] 370634.085 633117.096
[17,] 212022.909 370634.085
[18,] 257045.108 212022.909
[19,] -97859.522 257045.108
[20,] 713058.296 -97859.522
[21,] 164350.491 713058.296
[22,] 519750.260 164350.491
[23,] -97870.335 519750.260
[24,] -190459.271 -97870.335
[25,] -258595.195 -190459.271
[26,] 415969.927 -258595.195
[27,] -608210.939 415969.927
[28,] -109831.837 -608210.939
[29,] -11252.797 -109831.837
[30,] -208968.725 -11252.797
[31,] -259515.312 -208968.725
[32,] -407440.176 -259515.312
[33,] -417992.914 -407440.176
[34,] -231982.145 -417992.914
[35,] -877749.447 -231982.145
[36,] -67614.913 -877749.447
[37,] -161722.783 -67614.913
[38,] 128145.511 -161722.783
[39,] -461260.720 128145.511
[40,] 250036.215 -461260.720
[41,] -434512.183 250036.215
[42,] -243868.920 -434512.183
[43,] -259364.892 -243868.920
[44,] -61265.061 -259364.892
[45,] -235014.305 -61265.061
[46,] -706756.316 -235014.305
[47,] 385455.358 -706756.316
[48,] 338256.940 385455.358
[49,] 109289.497 338256.940
[50,] 38904.419 109289.497
[51,] 704602.269 38904.419
[52,] -186424.578 704602.269
[53,] 6323.156 -186424.578
[54,] 485428.480 6323.156
[55,] 464781.987 485428.480
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -337410.474 -683405.991
2 -747747.286 -337410.474
3 -268247.707 -747747.286
4 -324413.885 -268247.707
5 227418.914 -324413.885
6 -289635.943 227418.914
7 151957.740 -289635.943
8 -244353.058 151957.740
9 488656.727 -244353.058
10 418988.201 488656.727
11 590164.425 418988.201
12 603223.234 590164.425
13 648438.955 603223.234
14 164727.429 648438.955
15 633117.096 164727.429
16 370634.085 633117.096
17 212022.909 370634.085
18 257045.108 212022.909
19 -97859.522 257045.108
20 713058.296 -97859.522
21 164350.491 713058.296
22 519750.260 164350.491
23 -97870.335 519750.260
24 -190459.271 -97870.335
25 -258595.195 -190459.271
26 415969.927 -258595.195
27 -608210.939 415969.927
28 -109831.837 -608210.939
29 -11252.797 -109831.837
30 -208968.725 -11252.797
31 -259515.312 -208968.725
32 -407440.176 -259515.312
33 -417992.914 -407440.176
34 -231982.145 -417992.914
35 -877749.447 -231982.145
36 -67614.913 -877749.447
37 -161722.783 -67614.913
38 128145.511 -161722.783
39 -461260.720 128145.511
40 250036.215 -461260.720
41 -434512.183 250036.215
42 -243868.920 -434512.183
43 -259364.892 -243868.920
44 -61265.061 -259364.892
45 -235014.305 -61265.061
46 -706756.316 -235014.305
47 385455.358 -706756.316
48 338256.940 385455.358
49 109289.497 338256.940
50 38904.419 109289.497
51 704602.269 38904.419
52 -186424.578 704602.269
53 6323.156 -186424.578
54 485428.480 6323.156
55 464781.987 485428.480
> 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/7agay1290778168.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/8agay1290778168.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/9agay1290778168.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/10dhce1290778169.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/11zhs21290778169.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/12209q1290778169.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/13ysoz1290778169.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/14rjok1290778169.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/15u1m71290778169.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/168bkg1290778169.tab")
+ }
>
> try(system("convert tmp/1eouq1290778168.ps tmp/1eouq1290778168.png",intern=TRUE))
character(0)
> try(system("convert tmp/2eouq1290778168.ps tmp/2eouq1290778168.png",intern=TRUE))
character(0)
> try(system("convert tmp/3pxta1290778168.ps tmp/3pxta1290778168.png",intern=TRUE))
character(0)
> try(system("convert tmp/4pxta1290778168.ps tmp/4pxta1290778168.png",intern=TRUE))
character(0)
> try(system("convert tmp/5pxta1290778168.ps tmp/5pxta1290778168.png",intern=TRUE))
character(0)
> try(system("convert tmp/6iotd1290778168.ps tmp/6iotd1290778168.png",intern=TRUE))
character(0)
> try(system("convert tmp/7agay1290778168.ps tmp/7agay1290778168.png",intern=TRUE))
character(0)
> try(system("convert tmp/8agay1290778168.ps tmp/8agay1290778168.png",intern=TRUE))
character(0)
> try(system("convert tmp/9agay1290778168.ps tmp/9agay1290778168.png",intern=TRUE))
character(0)
> try(system("convert tmp/10dhce1290778169.ps tmp/10dhce1290778169.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
2.315 1.576 7.154