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(8823
+ ,0
+ ,9051
+ ,8776
+ ,0
+ ,8823
+ ,8255
+ ,0
+ ,8776
+ ,7969
+ ,0
+ ,8255
+ ,8758
+ ,0
+ ,7969
+ ,8693
+ ,0
+ ,8758
+ ,8271
+ ,0
+ ,8693
+ ,7790
+ ,0
+ ,8271
+ ,7769
+ ,0
+ ,7790
+ ,8170
+ ,0
+ ,7769
+ ,8209
+ ,0
+ ,8170
+ ,9395
+ ,0
+ ,8209
+ ,9260
+ ,0
+ ,9395
+ ,9018
+ ,0
+ ,9260
+ ,8501
+ ,0
+ ,9018
+ ,8500
+ ,0
+ ,8501
+ ,9649
+ ,0
+ ,8500
+ ,9319
+ ,0
+ ,9649
+ ,8830
+ ,0
+ ,9319
+ ,8436
+ ,0
+ ,8830
+ ,8169
+ ,0
+ ,8436
+ ,8269
+ ,0
+ ,8169
+ ,7945
+ ,0
+ ,8269
+ ,9144
+ ,0
+ ,7945
+ ,8770
+ ,0
+ ,9144
+ ,8834
+ ,0
+ ,8770
+ ,7837
+ ,0
+ ,8834
+ ,7792
+ ,0
+ ,7837
+ ,8616
+ ,0
+ ,7792
+ ,8518
+ ,0
+ ,8616
+ ,7940
+ ,0
+ ,8518
+ ,7545
+ ,0
+ ,7940
+ ,7531
+ ,0
+ ,7545
+ ,7665
+ ,0
+ ,7531
+ ,7599
+ ,0
+ ,7665
+ ,8444
+ ,0
+ ,7599
+ ,8549
+ ,0
+ ,8444
+ ,7986
+ ,0
+ ,8549
+ ,7335
+ ,0
+ ,7986
+ ,7287
+ ,0
+ ,7335
+ ,7870
+ ,0
+ ,7287
+ ,7839
+ ,0
+ ,7870
+ ,7327
+ ,0
+ ,7839
+ ,7259
+ ,0
+ ,7327
+ ,6964
+ ,0
+ ,7259
+ ,7271
+ ,0
+ ,6964
+ ,6956
+ ,0
+ ,7271
+ ,7608
+ ,0
+ ,6956
+ ,7692
+ ,0
+ ,7608
+ ,7255
+ ,0
+ ,7692
+ ,6804
+ ,0
+ ,7255
+ ,6655
+ ,0
+ ,6804
+ ,7341
+ ,0
+ ,6655
+ ,7602
+ ,0
+ ,7341
+ ,7086
+ ,0
+ ,7602
+ ,6625
+ ,0
+ ,7086
+ ,6272
+ ,0
+ ,6625
+ ,6576
+ ,0
+ ,6272
+ ,6491
+ ,0
+ ,6576
+ ,7649
+ ,0
+ ,6491
+ ,7400
+ ,0
+ ,7649
+ ,6913
+ ,0
+ ,7400
+ ,6532
+ ,0
+ ,6913
+ ,6486
+ ,0
+ ,6532
+ ,7295
+ ,0
+ ,6486
+ ,7556
+ ,0
+ ,7295
+ ,7088
+ ,1
+ ,7556
+ ,6952
+ ,1
+ ,7088
+ ,6773
+ ,1
+ ,6952
+ ,6917
+ ,1
+ ,6773
+ ,7371
+ ,1
+ ,6917
+ ,8221
+ ,1
+ ,7371
+ ,7953
+ ,1
+ ,8221
+ ,8027
+ ,1
+ ,7953
+ ,7287
+ ,1
+ ,8027
+ ,8076
+ ,1
+ ,7287
+ ,8933
+ ,1
+ ,8076
+ ,9433
+ ,1
+ ,8933
+ ,9479
+ ,1
+ ,9433
+ ,9199
+ ,1
+ ,9479
+ ,9469
+ ,1
+ ,9199
+ ,10015
+ ,1
+ ,9469
+ ,10999
+ ,1
+ ,10015
+ ,13009
+ ,1
+ ,10999
+ ,13699
+ ,1
+ ,13009
+ ,13895
+ ,1
+ ,13699
+ ,13248
+ ,1
+ ,13895
+ ,13973
+ ,1
+ ,13248
+ ,15095
+ ,1
+ ,13973
+ ,15201
+ ,1
+ ,15095
+ ,14823
+ ,1
+ ,15201
+ ,14538
+ ,1
+ ,14823
+ ,14547
+ ,1
+ ,14538
+ ,14407
+ ,1
+ ,14547)
+ ,dim=c(3
+ ,94)
+ ,dimnames=list(c('Y'
+ ,'X'
+ ,'Y1')
+ ,1:94))
> y <- array(NA,dim=c(3,94),dimnames=list(c('Y','X','Y1'),1:94))
> 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 M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 8823 0 9051 1 0 0 0 0 0 0 0 0 0 0 1
2 8776 0 8823 0 1 0 0 0 0 0 0 0 0 0 2
3 8255 0 8776 0 0 1 0 0 0 0 0 0 0 0 3
4 7969 0 8255 0 0 0 1 0 0 0 0 0 0 0 4
5 8758 0 7969 0 0 0 0 1 0 0 0 0 0 0 5
6 8693 0 8758 0 0 0 0 0 1 0 0 0 0 0 6
7 8271 0 8693 0 0 0 0 0 0 1 0 0 0 0 7
8 7790 0 8271 0 0 0 0 0 0 0 1 0 0 0 8
9 7769 0 7790 0 0 0 0 0 0 0 0 1 0 0 9
10 8170 0 7769 0 0 0 0 0 0 0 0 0 1 0 10
11 8209 0 8170 0 0 0 0 0 0 0 0 0 0 1 11
12 9395 0 8209 0 0 0 0 0 0 0 0 0 0 0 12
13 9260 0 9395 1 0 0 0 0 0 0 0 0 0 0 13
14 9018 0 9260 0 1 0 0 0 0 0 0 0 0 0 14
15 8501 0 9018 0 0 1 0 0 0 0 0 0 0 0 15
16 8500 0 8501 0 0 0 1 0 0 0 0 0 0 0 16
17 9649 0 8500 0 0 0 0 1 0 0 0 0 0 0 17
18 9319 0 9649 0 0 0 0 0 1 0 0 0 0 0 18
19 8830 0 9319 0 0 0 0 0 0 1 0 0 0 0 19
20 8436 0 8830 0 0 0 0 0 0 0 1 0 0 0 20
21 8169 0 8436 0 0 0 0 0 0 0 0 1 0 0 21
22 8269 0 8169 0 0 0 0 0 0 0 0 0 1 0 22
23 7945 0 8269 0 0 0 0 0 0 0 0 0 0 1 23
24 9144 0 7945 0 0 0 0 0 0 0 0 0 0 0 24
25 8770 0 9144 1 0 0 0 0 0 0 0 0 0 0 25
26 8834 0 8770 0 1 0 0 0 0 0 0 0 0 0 26
27 7837 0 8834 0 0 1 0 0 0 0 0 0 0 0 27
28 7792 0 7837 0 0 0 1 0 0 0 0 0 0 0 28
29 8616 0 7792 0 0 0 0 1 0 0 0 0 0 0 29
30 8518 0 8616 0 0 0 0 0 1 0 0 0 0 0 30
31 7940 0 8518 0 0 0 0 0 0 1 0 0 0 0 31
32 7545 0 7940 0 0 0 0 0 0 0 1 0 0 0 32
33 7531 0 7545 0 0 0 0 0 0 0 0 1 0 0 33
34 7665 0 7531 0 0 0 0 0 0 0 0 0 1 0 34
35 7599 0 7665 0 0 0 0 0 0 0 0 0 0 1 35
36 8444 0 7599 0 0 0 0 0 0 0 0 0 0 0 36
37 8549 0 8444 1 0 0 0 0 0 0 0 0 0 0 37
38 7986 0 8549 0 1 0 0 0 0 0 0 0 0 0 38
39 7335 0 7986 0 0 1 0 0 0 0 0 0 0 0 39
40 7287 0 7335 0 0 0 1 0 0 0 0 0 0 0 40
41 7870 0 7287 0 0 0 0 1 0 0 0 0 0 0 41
42 7839 0 7870 0 0 0 0 0 1 0 0 0 0 0 42
43 7327 0 7839 0 0 0 0 0 0 1 0 0 0 0 43
44 7259 0 7327 0 0 0 0 0 0 0 1 0 0 0 44
45 6964 0 7259 0 0 0 0 0 0 0 0 1 0 0 45
46 7271 0 6964 0 0 0 0 0 0 0 0 0 1 0 46
47 6956 0 7271 0 0 0 0 0 0 0 0 0 0 1 47
48 7608 0 6956 0 0 0 0 0 0 0 0 0 0 0 48
49 7692 0 7608 1 0 0 0 0 0 0 0 0 0 0 49
50 7255 0 7692 0 1 0 0 0 0 0 0 0 0 0 50
51 6804 0 7255 0 0 1 0 0 0 0 0 0 0 0 51
52 6655 0 6804 0 0 0 1 0 0 0 0 0 0 0 52
53 7341 0 6655 0 0 0 0 1 0 0 0 0 0 0 53
54 7602 0 7341 0 0 0 0 0 1 0 0 0 0 0 54
55 7086 0 7602 0 0 0 0 0 0 1 0 0 0 0 55
56 6625 0 7086 0 0 0 0 0 0 0 1 0 0 0 56
57 6272 0 6625 0 0 0 0 0 0 0 0 1 0 0 57
58 6576 0 6272 0 0 0 0 0 0 0 0 0 1 0 58
59 6491 0 6576 0 0 0 0 0 0 0 0 0 0 1 59
60 7649 0 6491 0 0 0 0 0 0 0 0 0 0 0 60
61 7400 0 7649 1 0 0 0 0 0 0 0 0 0 0 61
62 6913 0 7400 0 1 0 0 0 0 0 0 0 0 0 62
63 6532 0 6913 0 0 1 0 0 0 0 0 0 0 0 63
64 6486 0 6532 0 0 0 1 0 0 0 0 0 0 0 64
65 7295 0 6486 0 0 0 0 1 0 0 0 0 0 0 65
66 7556 0 7295 0 0 0 0 0 1 0 0 0 0 0 66
67 7088 1 7556 0 0 0 0 0 0 1 0 0 0 0 67
68 6952 1 7088 0 0 0 0 0 0 0 1 0 0 0 68
69 6773 1 6952 0 0 0 0 0 0 0 0 1 0 0 69
70 6917 1 6773 0 0 0 0 0 0 0 0 0 1 0 70
71 7371 1 6917 0 0 0 0 0 0 0 0 0 0 1 71
72 8221 1 7371 0 0 0 0 0 0 0 0 0 0 0 72
73 7953 1 8221 1 0 0 0 0 0 0 0 0 0 0 73
74 8027 1 7953 0 1 0 0 0 0 0 0 0 0 0 74
75 7287 1 8027 0 0 1 0 0 0 0 0 0 0 0 75
76 8076 1 7287 0 0 0 1 0 0 0 0 0 0 0 76
77 8933 1 8076 0 0 0 0 1 0 0 0 0 0 0 77
78 9433 1 8933 0 0 0 0 0 1 0 0 0 0 0 78
79 9479 1 9433 0 0 0 0 0 0 1 0 0 0 0 79
80 9199 1 9479 0 0 0 0 0 0 0 1 0 0 0 80
81 9469 1 9199 0 0 0 0 0 0 0 0 1 0 0 81
82 10015 1 9469 0 0 0 0 0 0 0 0 0 1 0 82
83 10999 1 10015 0 0 0 0 0 0 0 0 0 0 1 83
84 13009 1 10999 0 0 0 0 0 0 0 0 0 0 0 84
85 13699 1 13009 1 0 0 0 0 0 0 0 0 0 0 85
86 13895 1 13699 0 1 0 0 0 0 0 0 0 0 0 86
87 13248 1 13895 0 0 1 0 0 0 0 0 0 0 0 87
88 13973 1 13248 0 0 0 1 0 0 0 0 0 0 0 88
89 15095 1 13973 0 0 0 0 1 0 0 0 0 0 0 89
90 15201 1 15095 0 0 0 0 0 1 0 0 0 0 0 90
91 14823 1 15201 0 0 0 0 0 0 1 0 0 0 0 91
92 14538 1 14823 0 0 0 0 0 0 0 1 0 0 0 92
93 14547 1 14538 0 0 0 0 0 0 0 0 1 0 0 93
94 14407 1 14547 0 0 0 0 0 0 0 0 0 1 0 94
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X Y1 M1 M2 M3
1.029e+03 2.936e+02 1.002e+00 -1.166e+03 -1.300e+03 -1.732e+03
M4 M5 M6 M7 M8 M9
-1.001e+03 -2.663e+02 -1.045e+03 -1.572e+03 -1.469e+03 -1.262e+03
M10 M11 t
-9.313e+02 -1.030e+03 5.568e-02
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-560.173 -165.648 -8.735 170.178 670.624
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.029e+03 1.631e+02 6.311 1.50e-08 ***
X 2.936e+02 1.065e+02 2.756 0.00726 **
Y1 1.002e+00 1.567e-02 63.921 < 2e-16 ***
M1 -1.166e+03 1.361e+02 -8.568 6.78e-13 ***
M2 -1.300e+03 1.360e+02 -9.558 7.95e-15 ***
M3 -1.732e+03 1.356e+02 -12.774 < 2e-16 ***
M4 -1.001e+03 1.349e+02 -7.423 1.15e-10 ***
M5 -2.663e+02 1.350e+02 -1.973 0.05197 .
M6 -1.045e+03 1.364e+02 -7.657 4.05e-11 ***
M7 -1.572e+03 1.360e+02 -11.557 < 2e-16 ***
M8 -1.469e+03 1.353e+02 -10.858 < 2e-16 ***
M9 -1.262e+03 1.350e+02 -9.353 2.00e-14 ***
M10 -9.313e+02 1.349e+02 -6.905 1.14e-09 ***
M11 -1.030e+03 1.391e+02 -7.404 1.25e-10 ***
t 5.568e-02 1.630e+00 0.034 0.97284
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 260.3 on 79 degrees of freedom
Multiple R-squared: 0.9883, Adjusted R-squared: 0.9862
F-statistic: 474.7 on 14 and 79 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.3345797563 0.6691595127 0.6654202
[2,] 0.1854594410 0.3709188820 0.8145406
[3,] 0.1005984568 0.2011969137 0.8994015
[4,] 0.0647015880 0.1294031760 0.9352984
[5,] 0.0778680345 0.1557360691 0.9221320
[6,] 0.1226563693 0.2453127387 0.8773436
[7,] 0.0755048143 0.1510096286 0.9244952
[8,] 0.0511860144 0.1023720287 0.9488140
[9,] 0.0504662972 0.1009325943 0.9495337
[10,] 0.0918546752 0.1837093504 0.9081453
[11,] 0.0649383860 0.1298767719 0.9350616
[12,] 0.0437968419 0.0875936839 0.9562032
[13,] 0.0290814662 0.0581629323 0.9709185
[14,] 0.0178264904 0.0356529808 0.9821735
[15,] 0.0104107296 0.0208214591 0.9895893
[16,] 0.0077195687 0.0154391374 0.9922804
[17,] 0.0047254712 0.0094509424 0.9952745
[18,] 0.0026464351 0.0052928703 0.9973536
[19,] 0.0031707188 0.0063414377 0.9968293
[20,] 0.0053661541 0.0107323082 0.9946338
[21,] 0.0107348834 0.0214697668 0.9892651
[22,] 0.0066290222 0.0132580445 0.9933710
[23,] 0.0039450798 0.0078901597 0.9960549
[24,] 0.0040371316 0.0080742633 0.9959629
[25,] 0.0025244748 0.0050489496 0.9974755
[26,] 0.0014121251 0.0028242502 0.9985879
[27,] 0.0027883247 0.0055766495 0.9972117
[28,] 0.0019380371 0.0038760741 0.9980620
[29,] 0.0020410153 0.0040820305 0.9979590
[30,] 0.0017452768 0.0034905536 0.9982547
[31,] 0.0032891132 0.0065782263 0.9967109
[32,] 0.0032158502 0.0064317004 0.9967841
[33,] 0.0024029967 0.0048059935 0.9975970
[34,] 0.0031033648 0.0062067296 0.9968966
[35,] 0.0020314691 0.0040629381 0.9979685
[36,] 0.0013192620 0.0026385241 0.9986807
[37,] 0.0021719729 0.0043439458 0.9978280
[38,] 0.0013724842 0.0027449684 0.9986275
[39,] 0.0009613348 0.0019226696 0.9990387
[40,] 0.0006546190 0.0013092379 0.9993454
[41,] 0.0014776683 0.0029553365 0.9985223
[42,] 0.0009851292 0.0019702584 0.9990149
[43,] 0.0006757220 0.0013514439 0.9993243
[44,] 0.0003652265 0.0007304531 0.9996348
[45,] 0.0003160717 0.0006321434 0.9996839
[46,] 0.0007547809 0.0015095618 0.9992452
[47,] 0.0009465841 0.0018931683 0.9990534
[48,] 0.0004946894 0.0009893788 0.9995053
[49,] 0.0004239475 0.0008478950 0.9995761
[50,] 0.0002010707 0.0004021414 0.9997989
[51,] 0.0003532327 0.0007064653 0.9996468
[52,] 0.0002073429 0.0004146858 0.9997927
[53,] 0.0013731566 0.0027463132 0.9986268
[54,] 0.0043793792 0.0087587583 0.9956206
[55,] 0.0034451057 0.0068902114 0.9965549
[56,] 0.1730886094 0.3461772187 0.8269114
[57,] 0.1661799187 0.3323598373 0.8338201
[58,] 0.1290443850 0.2580887699 0.8709556
[59,] 0.1323557905 0.2647115811 0.8676442
> postscript(file="/var/www/html/rcomp/tmp/1jc2b1260040089.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/2mr7h1260040089.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/305hp1260040089.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/427uc1260040089.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/5jebp1260040089.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 = 94
Frequency = 1
1 2 3 4 5 6
-105.356356 209.309611 167.970553 -327.677441 12.970557 -64.052635
7 8 9 10 11 12
106.001862 -55.111290 198.909395 290.022446 26.226404 142.892842
13 14 15 16 17 18
-13.577077 12.939491 170.913688 -43.740731 371.449432 -331.152079
19 20 21 22 23 24
37.328125 30.322604 -48.796471 -12.288236 -337.600739 155.648837
25 26 27 28 29 30
-252.841966 319.058541 -309.458832 -87.342160 46.918686 -98.160731
31 32 33 34 35 36
-51.053221 30.084231 204.966762 22.068568 -79.298554 -198.463440
37 38 39 40 41 42
226.614437 -308.254541 37.235324 -90.203832 -193.938167 -30.630432
43 44 45 46 47 48
15.369446 357.400874 -76.241901 195.311314 -328.333722 -391.098604
49 50 51 52 53 54
206.289316 -181.545927 237.741527 -191.018917 -90.591001 261.551270
55 56 57 58 59 60
11.082067 -35.880080 -133.891523 192.754866 -97.885351 114.980285
61 62 63 64 65 66
-127.444653 -231.744951 307.622824 -88.250071 32.012382 260.957061
67 68 69 70 71 72
-235.142478 -5.181734 -254.715271 -262.348402 146.268411 -488.731824
73 74 75 76 77 78
-441.661979 34.068245 -347.465192 451.238625 -216.840301 203.027269
79 80 81 82 83 84
275.174127 -153.690840 190.006951 134.652528 670.623552 664.771903
85 86 87 88 89 90
507.978277 146.169531 -264.559893 376.994526 38.018412 -201.539725
91 92 93 94
-158.759929 -167.943765 -80.237942 -560.173084
> postscript(file="/var/www/html/rcomp/tmp/6z0u41260040089.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 = 94
Frequency = 1
lag(myerror, k = 1) myerror
0 -105.356356 NA
1 209.309611 -105.356356
2 167.970553 209.309611
3 -327.677441 167.970553
4 12.970557 -327.677441
5 -64.052635 12.970557
6 106.001862 -64.052635
7 -55.111290 106.001862
8 198.909395 -55.111290
9 290.022446 198.909395
10 26.226404 290.022446
11 142.892842 26.226404
12 -13.577077 142.892842
13 12.939491 -13.577077
14 170.913688 12.939491
15 -43.740731 170.913688
16 371.449432 -43.740731
17 -331.152079 371.449432
18 37.328125 -331.152079
19 30.322604 37.328125
20 -48.796471 30.322604
21 -12.288236 -48.796471
22 -337.600739 -12.288236
23 155.648837 -337.600739
24 -252.841966 155.648837
25 319.058541 -252.841966
26 -309.458832 319.058541
27 -87.342160 -309.458832
28 46.918686 -87.342160
29 -98.160731 46.918686
30 -51.053221 -98.160731
31 30.084231 -51.053221
32 204.966762 30.084231
33 22.068568 204.966762
34 -79.298554 22.068568
35 -198.463440 -79.298554
36 226.614437 -198.463440
37 -308.254541 226.614437
38 37.235324 -308.254541
39 -90.203832 37.235324
40 -193.938167 -90.203832
41 -30.630432 -193.938167
42 15.369446 -30.630432
43 357.400874 15.369446
44 -76.241901 357.400874
45 195.311314 -76.241901
46 -328.333722 195.311314
47 -391.098604 -328.333722
48 206.289316 -391.098604
49 -181.545927 206.289316
50 237.741527 -181.545927
51 -191.018917 237.741527
52 -90.591001 -191.018917
53 261.551270 -90.591001
54 11.082067 261.551270
55 -35.880080 11.082067
56 -133.891523 -35.880080
57 192.754866 -133.891523
58 -97.885351 192.754866
59 114.980285 -97.885351
60 -127.444653 114.980285
61 -231.744951 -127.444653
62 307.622824 -231.744951
63 -88.250071 307.622824
64 32.012382 -88.250071
65 260.957061 32.012382
66 -235.142478 260.957061
67 -5.181734 -235.142478
68 -254.715271 -5.181734
69 -262.348402 -254.715271
70 146.268411 -262.348402
71 -488.731824 146.268411
72 -441.661979 -488.731824
73 34.068245 -441.661979
74 -347.465192 34.068245
75 451.238625 -347.465192
76 -216.840301 451.238625
77 203.027269 -216.840301
78 275.174127 203.027269
79 -153.690840 275.174127
80 190.006951 -153.690840
81 134.652528 190.006951
82 670.623552 134.652528
83 664.771903 670.623552
84 507.978277 664.771903
85 146.169531 507.978277
86 -264.559893 146.169531
87 376.994526 -264.559893
88 38.018412 376.994526
89 -201.539725 38.018412
90 -158.759929 -201.539725
91 -167.943765 -158.759929
92 -80.237942 -167.943765
93 -560.173084 -80.237942
94 NA -560.173084
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 209.309611 -105.356356
[2,] 167.970553 209.309611
[3,] -327.677441 167.970553
[4,] 12.970557 -327.677441
[5,] -64.052635 12.970557
[6,] 106.001862 -64.052635
[7,] -55.111290 106.001862
[8,] 198.909395 -55.111290
[9,] 290.022446 198.909395
[10,] 26.226404 290.022446
[11,] 142.892842 26.226404
[12,] -13.577077 142.892842
[13,] 12.939491 -13.577077
[14,] 170.913688 12.939491
[15,] -43.740731 170.913688
[16,] 371.449432 -43.740731
[17,] -331.152079 371.449432
[18,] 37.328125 -331.152079
[19,] 30.322604 37.328125
[20,] -48.796471 30.322604
[21,] -12.288236 -48.796471
[22,] -337.600739 -12.288236
[23,] 155.648837 -337.600739
[24,] -252.841966 155.648837
[25,] 319.058541 -252.841966
[26,] -309.458832 319.058541
[27,] -87.342160 -309.458832
[28,] 46.918686 -87.342160
[29,] -98.160731 46.918686
[30,] -51.053221 -98.160731
[31,] 30.084231 -51.053221
[32,] 204.966762 30.084231
[33,] 22.068568 204.966762
[34,] -79.298554 22.068568
[35,] -198.463440 -79.298554
[36,] 226.614437 -198.463440
[37,] -308.254541 226.614437
[38,] 37.235324 -308.254541
[39,] -90.203832 37.235324
[40,] -193.938167 -90.203832
[41,] -30.630432 -193.938167
[42,] 15.369446 -30.630432
[43,] 357.400874 15.369446
[44,] -76.241901 357.400874
[45,] 195.311314 -76.241901
[46,] -328.333722 195.311314
[47,] -391.098604 -328.333722
[48,] 206.289316 -391.098604
[49,] -181.545927 206.289316
[50,] 237.741527 -181.545927
[51,] -191.018917 237.741527
[52,] -90.591001 -191.018917
[53,] 261.551270 -90.591001
[54,] 11.082067 261.551270
[55,] -35.880080 11.082067
[56,] -133.891523 -35.880080
[57,] 192.754866 -133.891523
[58,] -97.885351 192.754866
[59,] 114.980285 -97.885351
[60,] -127.444653 114.980285
[61,] -231.744951 -127.444653
[62,] 307.622824 -231.744951
[63,] -88.250071 307.622824
[64,] 32.012382 -88.250071
[65,] 260.957061 32.012382
[66,] -235.142478 260.957061
[67,] -5.181734 -235.142478
[68,] -254.715271 -5.181734
[69,] -262.348402 -254.715271
[70,] 146.268411 -262.348402
[71,] -488.731824 146.268411
[72,] -441.661979 -488.731824
[73,] 34.068245 -441.661979
[74,] -347.465192 34.068245
[75,] 451.238625 -347.465192
[76,] -216.840301 451.238625
[77,] 203.027269 -216.840301
[78,] 275.174127 203.027269
[79,] -153.690840 275.174127
[80,] 190.006951 -153.690840
[81,] 134.652528 190.006951
[82,] 670.623552 134.652528
[83,] 664.771903 670.623552
[84,] 507.978277 664.771903
[85,] 146.169531 507.978277
[86,] -264.559893 146.169531
[87,] 376.994526 -264.559893
[88,] 38.018412 376.994526
[89,] -201.539725 38.018412
[90,] -158.759929 -201.539725
[91,] -167.943765 -158.759929
[92,] -80.237942 -167.943765
[93,] -560.173084 -80.237942
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 209.309611 -105.356356
2 167.970553 209.309611
3 -327.677441 167.970553
4 12.970557 -327.677441
5 -64.052635 12.970557
6 106.001862 -64.052635
7 -55.111290 106.001862
8 198.909395 -55.111290
9 290.022446 198.909395
10 26.226404 290.022446
11 142.892842 26.226404
12 -13.577077 142.892842
13 12.939491 -13.577077
14 170.913688 12.939491
15 -43.740731 170.913688
16 371.449432 -43.740731
17 -331.152079 371.449432
18 37.328125 -331.152079
19 30.322604 37.328125
20 -48.796471 30.322604
21 -12.288236 -48.796471
22 -337.600739 -12.288236
23 155.648837 -337.600739
24 -252.841966 155.648837
25 319.058541 -252.841966
26 -309.458832 319.058541
27 -87.342160 -309.458832
28 46.918686 -87.342160
29 -98.160731 46.918686
30 -51.053221 -98.160731
31 30.084231 -51.053221
32 204.966762 30.084231
33 22.068568 204.966762
34 -79.298554 22.068568
35 -198.463440 -79.298554
36 226.614437 -198.463440
37 -308.254541 226.614437
38 37.235324 -308.254541
39 -90.203832 37.235324
40 -193.938167 -90.203832
41 -30.630432 -193.938167
42 15.369446 -30.630432
43 357.400874 15.369446
44 -76.241901 357.400874
45 195.311314 -76.241901
46 -328.333722 195.311314
47 -391.098604 -328.333722
48 206.289316 -391.098604
49 -181.545927 206.289316
50 237.741527 -181.545927
51 -191.018917 237.741527
52 -90.591001 -191.018917
53 261.551270 -90.591001
54 11.082067 261.551270
55 -35.880080 11.082067
56 -133.891523 -35.880080
57 192.754866 -133.891523
58 -97.885351 192.754866
59 114.980285 -97.885351
60 -127.444653 114.980285
61 -231.744951 -127.444653
62 307.622824 -231.744951
63 -88.250071 307.622824
64 32.012382 -88.250071
65 260.957061 32.012382
66 -235.142478 260.957061
67 -5.181734 -235.142478
68 -254.715271 -5.181734
69 -262.348402 -254.715271
70 146.268411 -262.348402
71 -488.731824 146.268411
72 -441.661979 -488.731824
73 34.068245 -441.661979
74 -347.465192 34.068245
75 451.238625 -347.465192
76 -216.840301 451.238625
77 203.027269 -216.840301
78 275.174127 203.027269
79 -153.690840 275.174127
80 190.006951 -153.690840
81 134.652528 190.006951
82 670.623552 134.652528
83 664.771903 670.623552
84 507.978277 664.771903
85 146.169531 507.978277
86 -264.559893 146.169531
87 376.994526 -264.559893
88 38.018412 376.994526
89 -201.539725 38.018412
90 -158.759929 -201.539725
91 -167.943765 -158.759929
92 -80.237942 -167.943765
93 -560.173084 -80.237942
> 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/7zmwj1260040089.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/8yyu81260040089.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/9h8ze1260040089.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/10prl21260040089.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/11xpn01260040089.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/12huss1260040089.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/130whs1260040090.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/14fr7e1260040090.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/15sn641260040090.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/16w7371260040090.tab")
+ }
>
> system("convert tmp/1jc2b1260040089.ps tmp/1jc2b1260040089.png")
> system("convert tmp/2mr7h1260040089.ps tmp/2mr7h1260040089.png")
> system("convert tmp/305hp1260040089.ps tmp/305hp1260040089.png")
> system("convert tmp/427uc1260040089.ps tmp/427uc1260040089.png")
> system("convert tmp/5jebp1260040089.ps tmp/5jebp1260040089.png")
> system("convert tmp/6z0u41260040089.ps tmp/6z0u41260040089.png")
> system("convert tmp/7zmwj1260040089.ps tmp/7zmwj1260040089.png")
> system("convert tmp/8yyu81260040089.ps tmp/8yyu81260040089.png")
> system("convert tmp/9h8ze1260040089.ps tmp/9h8ze1260040089.png")
> system("convert tmp/10prl21260040089.ps tmp/10prl21260040089.png")
>
>
> proc.time()
user system elapsed
2.913 1.611 3.392