R version 2.15.2 (2012-10-26) -- "Trick or Treat"
Copyright (C) 2012 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: i686-pc-linux-gnu (32-bit)
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(68.897,38.683,44.720,39.525,45.315,50.380,40.600,36.279,42.438,38.064,31.879,11.379,70.249,39.253,47.060,41.697,38.708,49.267,39.018,32.228,40.870,39.383,34.571,12.066,70.938,34.077,45.409,40.809,37.013,44.953,37.848,32.745,43.412,34.931,33.008,8.620,68.906,39.556,50.669,36.432,40.891,48.428,36.222,33.425,39.401,37.967,34.801,12.657,69.116,41.519,51.321,38.529,41.547,52.073,38.401,40.898,40.439,41.888,37.898,8.771,68.184,50.530,47.221,41.756,45.633,48.138,39.486,39.341,41.117,41.629,29.722,7.054,56.676,34.870,35.117,30.169,30.936,35.699,33.228,27.733,33.666,35.429,27.438,8.170,63.410,38.040,45.389,37.353,37.024,50.957,37.994,36.454,46.080,43.373,37.395,10.963,76.058,50.179,57.452,47.568,50.050,50.856,41.992,39.284,44.521,43.832,41.153,17.100),dim=c(1,108),dimnames=list(c('maandelijkseverkoopcijfers'),1:108))
> y <- array(NA,dim=c(1,108),dimnames=list(c('maandelijkseverkoopcijfers'),1:108))
> 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'
> 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, 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
maandelijkseverkoopcijfers M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 68.897 1 0 0 0 0 0 0 0 0 0 0 1
2 38.683 0 1 0 0 0 0 0 0 0 0 0 2
3 44.720 0 0 1 0 0 0 0 0 0 0 0 3
4 39.525 0 0 0 1 0 0 0 0 0 0 0 4
5 45.315 0 0 0 0 1 0 0 0 0 0 0 5
6 50.380 0 0 0 0 0 1 0 0 0 0 0 6
7 40.600 0 0 0 0 0 0 1 0 0 0 0 7
8 36.279 0 0 0 0 0 0 0 1 0 0 0 8
9 42.438 0 0 0 0 0 0 0 0 1 0 0 9
10 38.064 0 0 0 0 0 0 0 0 0 1 0 10
11 31.879 0 0 0 0 0 0 0 0 0 0 1 11
12 11.379 0 0 0 0 0 0 0 0 0 0 0 12
13 70.249 1 0 0 0 0 0 0 0 0 0 0 13
14 39.253 0 1 0 0 0 0 0 0 0 0 0 14
15 47.060 0 0 1 0 0 0 0 0 0 0 0 15
16 41.697 0 0 0 1 0 0 0 0 0 0 0 16
17 38.708 0 0 0 0 1 0 0 0 0 0 0 17
18 49.267 0 0 0 0 0 1 0 0 0 0 0 18
19 39.018 0 0 0 0 0 0 1 0 0 0 0 19
20 32.228 0 0 0 0 0 0 0 1 0 0 0 20
21 40.870 0 0 0 0 0 0 0 0 1 0 0 21
22 39.383 0 0 0 0 0 0 0 0 0 1 0 22
23 34.571 0 0 0 0 0 0 0 0 0 0 1 23
24 12.066 0 0 0 0 0 0 0 0 0 0 0 24
25 70.938 1 0 0 0 0 0 0 0 0 0 0 25
26 34.077 0 1 0 0 0 0 0 0 0 0 0 26
27 45.409 0 0 1 0 0 0 0 0 0 0 0 27
28 40.809 0 0 0 1 0 0 0 0 0 0 0 28
29 37.013 0 0 0 0 1 0 0 0 0 0 0 29
30 44.953 0 0 0 0 0 1 0 0 0 0 0 30
31 37.848 0 0 0 0 0 0 1 0 0 0 0 31
32 32.745 0 0 0 0 0 0 0 1 0 0 0 32
33 43.412 0 0 0 0 0 0 0 0 1 0 0 33
34 34.931 0 0 0 0 0 0 0 0 0 1 0 34
35 33.008 0 0 0 0 0 0 0 0 0 0 1 35
36 8.620 0 0 0 0 0 0 0 0 0 0 0 36
37 68.906 1 0 0 0 0 0 0 0 0 0 0 37
38 39.556 0 1 0 0 0 0 0 0 0 0 0 38
39 50.669 0 0 1 0 0 0 0 0 0 0 0 39
40 36.432 0 0 0 1 0 0 0 0 0 0 0 40
41 40.891 0 0 0 0 1 0 0 0 0 0 0 41
42 48.428 0 0 0 0 0 1 0 0 0 0 0 42
43 36.222 0 0 0 0 0 0 1 0 0 0 0 43
44 33.425 0 0 0 0 0 0 0 1 0 0 0 44
45 39.401 0 0 0 0 0 0 0 0 1 0 0 45
46 37.967 0 0 0 0 0 0 0 0 0 1 0 46
47 34.801 0 0 0 0 0 0 0 0 0 0 1 47
48 12.657 0 0 0 0 0 0 0 0 0 0 0 48
49 69.116 1 0 0 0 0 0 0 0 0 0 0 49
50 41.519 0 1 0 0 0 0 0 0 0 0 0 50
51 51.321 0 0 1 0 0 0 0 0 0 0 0 51
52 38.529 0 0 0 1 0 0 0 0 0 0 0 52
53 41.547 0 0 0 0 1 0 0 0 0 0 0 53
54 52.073 0 0 0 0 0 1 0 0 0 0 0 54
55 38.401 0 0 0 0 0 0 1 0 0 0 0 55
56 40.898 0 0 0 0 0 0 0 1 0 0 0 56
57 40.439 0 0 0 0 0 0 0 0 1 0 0 57
58 41.888 0 0 0 0 0 0 0 0 0 1 0 58
59 37.898 0 0 0 0 0 0 0 0 0 0 1 59
60 8.771 0 0 0 0 0 0 0 0 0 0 0 60
61 68.184 1 0 0 0 0 0 0 0 0 0 0 61
62 50.530 0 1 0 0 0 0 0 0 0 0 0 62
63 47.221 0 0 1 0 0 0 0 0 0 0 0 63
64 41.756 0 0 0 1 0 0 0 0 0 0 0 64
65 45.633 0 0 0 0 1 0 0 0 0 0 0 65
66 48.138 0 0 0 0 0 1 0 0 0 0 0 66
67 39.486 0 0 0 0 0 0 1 0 0 0 0 67
68 39.341 0 0 0 0 0 0 0 1 0 0 0 68
69 41.117 0 0 0 0 0 0 0 0 1 0 0 69
70 41.629 0 0 0 0 0 0 0 0 0 1 0 70
71 29.722 0 0 0 0 0 0 0 0 0 0 1 71
72 7.054 0 0 0 0 0 0 0 0 0 0 0 72
73 56.676 1 0 0 0 0 0 0 0 0 0 0 73
74 34.870 0 1 0 0 0 0 0 0 0 0 0 74
75 35.117 0 0 1 0 0 0 0 0 0 0 0 75
76 30.169 0 0 0 1 0 0 0 0 0 0 0 76
77 30.936 0 0 0 0 1 0 0 0 0 0 0 77
78 35.699 0 0 0 0 0 1 0 0 0 0 0 78
79 33.228 0 0 0 0 0 0 1 0 0 0 0 79
80 27.733 0 0 0 0 0 0 0 1 0 0 0 80
81 33.666 0 0 0 0 0 0 0 0 1 0 0 81
82 35.429 0 0 0 0 0 0 0 0 0 1 0 82
83 27.438 0 0 0 0 0 0 0 0 0 0 1 83
84 8.170 0 0 0 0 0 0 0 0 0 0 0 84
85 63.410 1 0 0 0 0 0 0 0 0 0 0 85
86 38.040 0 1 0 0 0 0 0 0 0 0 0 86
87 45.389 0 0 1 0 0 0 0 0 0 0 0 87
88 37.353 0 0 0 1 0 0 0 0 0 0 0 88
89 37.024 0 0 0 0 1 0 0 0 0 0 0 89
90 50.957 0 0 0 0 0 1 0 0 0 0 0 90
91 37.994 0 0 0 0 0 0 1 0 0 0 0 91
92 36.454 0 0 0 0 0 0 0 1 0 0 0 92
93 46.080 0 0 0 0 0 0 0 0 1 0 0 93
94 43.373 0 0 0 0 0 0 0 0 0 1 0 94
95 37.395 0 0 0 0 0 0 0 0 0 0 1 95
96 10.963 0 0 0 0 0 0 0 0 0 0 0 96
97 76.058 1 0 0 0 0 0 0 0 0 0 0 97
98 50.179 0 1 0 0 0 0 0 0 0 0 0 98
99 57.452 0 0 1 0 0 0 0 0 0 0 0 99
100 47.568 0 0 0 1 0 0 0 0 0 0 0 100
101 50.050 0 0 0 0 1 0 0 0 0 0 0 101
102 50.856 0 0 0 0 0 1 0 0 0 0 0 102
103 41.992 0 0 0 0 0 0 1 0 0 0 0 103
104 39.284 0 0 0 0 0 0 0 1 0 0 0 104
105 44.521 0 0 0 0 0 0 0 0 1 0 0 105
106 43.832 0 0 0 0 0 0 0 0 0 1 0 106
107 41.153 0 0 0 0 0 0 0 0 0 0 1 107
108 17.100 0 0 0 0 0 0 0 0 0 0 0 108
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) M1 M2 M3 M4 M5
9.64971 57.49722 30.17583 36.56310 28.70915 30.16620
M6 M7 M8 M9 M10 M11
37.21825 27.64852 24.69657 30.62896 28.89412 23.47228
t
0.01839
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-12.604 -1.896 0.425 2.526 9.564
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 9.64971 1.75635 5.494 3.28e-07 ***
M1 57.49722 2.17377 26.450 < 2e-16 ***
M2 30.17583 2.17278 13.888 < 2e-16 ***
M3 36.56310 2.17189 16.835 < 2e-16 ***
M4 28.70915 2.17109 13.223 < 2e-16 ***
M5 30.16620 2.17039 13.899 < 2e-16 ***
M6 37.21825 2.16978 17.153 < 2e-16 ***
M7 27.64852 2.16926 12.746 < 2e-16 ***
M8 24.69657 2.16884 11.387 < 2e-16 ***
M9 30.62896 2.16851 14.124 < 2e-16 ***
M10 28.89412 2.16827 13.326 < 2e-16 ***
M11 23.47228 2.16813 10.826 < 2e-16 ***
t 0.01839 0.01428 1.288 0.201
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 4.599 on 95 degrees of freedom
Multiple R-squared: 0.8907, Adjusted R-squared: 0.8769
F-statistic: 64.5 on 12 and 95 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,] 2.263813e-03 4.527626e-03 0.99773619
[2,] 1.158490e-01 2.316979e-01 0.88415103
[3,] 5.349870e-02 1.069974e-01 0.94650130
[4,] 2.377377e-02 4.754754e-02 0.97622623
[5,] 1.579998e-02 3.159996e-02 0.98420002
[6,] 6.408846e-03 1.281769e-02 0.99359115
[7,] 3.075192e-03 6.150384e-03 0.99692481
[8,] 1.994548e-03 3.989097e-03 0.99800545
[9,] 8.108036e-04 1.621607e-03 0.99918920
[10,] 3.895043e-04 7.790085e-04 0.99961050
[11,] 5.997859e-04 1.199572e-03 0.99940021
[12,] 2.331603e-04 4.663205e-04 0.99976684
[13,] 9.741752e-05 1.948350e-04 0.99990258
[14,] 1.117620e-04 2.235240e-04 0.99988824
[15,] 9.215325e-05 1.843065e-04 0.99990785
[16,] 3.630211e-05 7.260423e-05 0.99996370
[17,] 1.365663e-05 2.731326e-05 0.99998634
[18,] 1.093444e-05 2.186888e-05 0.99998907
[19,] 6.217003e-06 1.243401e-05 0.99999378
[20,] 2.524419e-06 5.048838e-06 0.99999748
[21,] 1.086645e-06 2.173289e-06 0.99999891
[22,] 4.365715e-07 8.731429e-07 0.99999956
[23,] 5.279681e-07 1.055936e-06 0.99999947
[24,] 2.801414e-06 5.602827e-06 0.99999720
[25,] 1.986435e-06 3.972871e-06 0.99999801
[26,] 9.343404e-07 1.868681e-06 0.99999907
[27,] 4.148649e-07 8.297298e-07 0.99999959
[28,] 1.949186e-07 3.898373e-07 0.99999981
[29,] 7.543370e-08 1.508674e-07 0.99999992
[30,] 3.320604e-08 6.641208e-08 0.99999997
[31,] 1.420794e-08 2.841588e-08 0.99999999
[32,] 8.245905e-09 1.649181e-08 0.99999999
[33,] 5.670386e-09 1.134077e-08 0.99999999
[34,] 2.403379e-09 4.806758e-09 1.00000000
[35,] 3.226994e-09 6.453988e-09 1.00000000
[36,] 7.616033e-09 1.523207e-08 0.99999999
[37,] 2.950129e-09 5.900257e-09 1.00000000
[38,] 1.296143e-09 2.592285e-09 1.00000000
[39,] 2.306811e-09 4.613621e-09 1.00000000
[40,] 9.306579e-10 1.861316e-09 1.00000000
[41,] 1.156481e-08 2.312962e-08 0.99999999
[42,] 5.396229e-09 1.079246e-08 0.99999999
[43,] 5.885188e-09 1.177038e-08 0.99999999
[44,] 1.004499e-08 2.008999e-08 0.99999999
[45,] 6.670816e-09 1.334163e-08 0.99999999
[46,] 5.414759e-09 1.082952e-08 0.99999999
[47,] 2.551450e-06 5.102900e-06 0.99999745
[48,] 2.176488e-06 4.352977e-06 0.99999782
[49,] 2.995818e-06 5.991635e-06 0.99999700
[50,] 1.737168e-05 3.474335e-05 0.99998263
[51,] 3.612057e-05 7.224114e-05 0.99996388
[52,] 7.904746e-05 1.580949e-04 0.99992095
[53,] 1.085560e-03 2.171120e-03 0.99891444
[54,] 3.537339e-03 7.074677e-03 0.99646266
[55,] 2.904011e-02 5.808022e-02 0.97095989
[56,] 6.945220e-02 1.389044e-01 0.93054780
[57,] 1.831990e-01 3.663980e-01 0.81680102
[58,] 4.212265e-01 8.424530e-01 0.57877349
[59,] 4.099970e-01 8.199939e-01 0.59000304
[60,] 6.459173e-01 7.081653e-01 0.35408267
[61,] 6.723867e-01 6.552265e-01 0.32761325
[62,] 7.042307e-01 5.915385e-01 0.29576925
[63,] 8.154580e-01 3.690840e-01 0.18454200
[64,] 7.807881e-01 4.384238e-01 0.21921189
[65,] 7.384364e-01 5.231271e-01 0.26156357
[66,] 6.998127e-01 6.003746e-01 0.30018731
[67,] 6.234431e-01 7.531137e-01 0.37655687
[68,] 5.717956e-01 8.564088e-01 0.42820440
[69,] 5.143634e-01 9.712733e-01 0.48563665
[70,] 5.373876e-01 9.252249e-01 0.46261244
[71,] 5.709630e-01 8.580739e-01 0.42903695
[72,] 6.474074e-01 7.051852e-01 0.35259262
[73,] 6.855237e-01 6.289527e-01 0.31447634
[74,] 9.606281e-01 7.874389e-02 0.03937195
[75,] 9.407570e-01 1.184860e-01 0.05924300
[76,] 8.833078e-01 2.333844e-01 0.11669218
[77,] 7.646624e-01 4.706753e-01 0.23533763
> postscript(file="/var/fisher/rcomp/tmp/1110l1355321774.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/fisher/rcomp/tmp/24xq51355321774.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/fisher/rcomp/tmp/3c8v61355321774.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/fisher/rcomp/tmp/4s5vj1355321774.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/fisher/rcomp/tmp/5j15z1355321774.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 = 108
Frequency = 1
1 2 3 4 5 6
1.73167778 -1.17932222 -1.54798889 1.09256667 5.40712222 3.40167778
7 8 9 10 11 12
3.17301111 1.78556667 1.99378889 -0.66376667 -1.44532222 1.50856667
13 14 15 16 17 18
2.86295278 -0.83004722 0.57128611 3.04384167 -1.42060278 2.06795278
19 20 21 22 23 24
1.37028611 -2.48615833 0.20506389 0.43450833 1.02595278 1.97484167
25 26 27 28 29 30
3.33122778 -6.22677222 -1.30043889 1.93511667 -3.33632778 -2.46677222
31 32 33 34 35 36
-0.02043889 -2.18988333 2.52633889 -4.23821667 -0.75777222 -1.69188333
37 38 39 40 41 42
1.07850278 -0.96849722 3.73883611 -2.66260833 0.32094722 0.78750278
43 44 45 46 47 48
-1.86716389 -1.73060833 -1.70538611 -1.42294167 0.81450278 2.12439167
49 50 51 52 53 54
1.06777778 0.77377778 4.17011111 -0.78633333 0.75622222 4.21177778
55 56 57 58 59 60
0.09111111 5.52166667 -0.88811111 2.27733333 3.69077778 -1.98233333
61 62 63 64 65 66
-0.08494722 9.56405278 -0.15061389 2.21994167 4.62149722 0.05605278
67 68 69 70 71 72
0.95538611 3.74394167 -0.43083611 1.79760833 -4.70594722 -3.92005833
73 74 75 76 77 78
-11.81367222 -6.31667222 -12.47533889 -9.58778333 -10.29622778 -12.60367222
79 80 81 82 83 84
-5.52333889 -8.08478333 -8.10256111 -4.62311667 -7.21067222 -3.02478333
85 86 87 88 89 90
-5.30039722 -3.36739722 -2.42406389 -2.62450833 -4.42895278 2.43360278
91 92 93 94 95 96
-0.97806389 0.41549167 4.09071389 3.10015833 2.52560278 -0.45250833
97 98 99 100 101 102
7.12687778 8.55087778 9.41821111 7.36976667 8.37632222 2.11187778
103 104 105 106 107 108
2.79921111 3.02476667 2.31098889 3.33843333 6.06287778 5.46376667
> postscript(file="/var/fisher/rcomp/tmp/6in0b1355321774.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 = 108
Frequency = 1
lag(myerror, k = 1) myerror
0 1.73167778 NA
1 -1.17932222 1.73167778
2 -1.54798889 -1.17932222
3 1.09256667 -1.54798889
4 5.40712222 1.09256667
5 3.40167778 5.40712222
6 3.17301111 3.40167778
7 1.78556667 3.17301111
8 1.99378889 1.78556667
9 -0.66376667 1.99378889
10 -1.44532222 -0.66376667
11 1.50856667 -1.44532222
12 2.86295278 1.50856667
13 -0.83004722 2.86295278
14 0.57128611 -0.83004722
15 3.04384167 0.57128611
16 -1.42060278 3.04384167
17 2.06795278 -1.42060278
18 1.37028611 2.06795278
19 -2.48615833 1.37028611
20 0.20506389 -2.48615833
21 0.43450833 0.20506389
22 1.02595278 0.43450833
23 1.97484167 1.02595278
24 3.33122778 1.97484167
25 -6.22677222 3.33122778
26 -1.30043889 -6.22677222
27 1.93511667 -1.30043889
28 -3.33632778 1.93511667
29 -2.46677222 -3.33632778
30 -0.02043889 -2.46677222
31 -2.18988333 -0.02043889
32 2.52633889 -2.18988333
33 -4.23821667 2.52633889
34 -0.75777222 -4.23821667
35 -1.69188333 -0.75777222
36 1.07850278 -1.69188333
37 -0.96849722 1.07850278
38 3.73883611 -0.96849722
39 -2.66260833 3.73883611
40 0.32094722 -2.66260833
41 0.78750278 0.32094722
42 -1.86716389 0.78750278
43 -1.73060833 -1.86716389
44 -1.70538611 -1.73060833
45 -1.42294167 -1.70538611
46 0.81450278 -1.42294167
47 2.12439167 0.81450278
48 1.06777778 2.12439167
49 0.77377778 1.06777778
50 4.17011111 0.77377778
51 -0.78633333 4.17011111
52 0.75622222 -0.78633333
53 4.21177778 0.75622222
54 0.09111111 4.21177778
55 5.52166667 0.09111111
56 -0.88811111 5.52166667
57 2.27733333 -0.88811111
58 3.69077778 2.27733333
59 -1.98233333 3.69077778
60 -0.08494722 -1.98233333
61 9.56405278 -0.08494722
62 -0.15061389 9.56405278
63 2.21994167 -0.15061389
64 4.62149722 2.21994167
65 0.05605278 4.62149722
66 0.95538611 0.05605278
67 3.74394167 0.95538611
68 -0.43083611 3.74394167
69 1.79760833 -0.43083611
70 -4.70594722 1.79760833
71 -3.92005833 -4.70594722
72 -11.81367222 -3.92005833
73 -6.31667222 -11.81367222
74 -12.47533889 -6.31667222
75 -9.58778333 -12.47533889
76 -10.29622778 -9.58778333
77 -12.60367222 -10.29622778
78 -5.52333889 -12.60367222
79 -8.08478333 -5.52333889
80 -8.10256111 -8.08478333
81 -4.62311667 -8.10256111
82 -7.21067222 -4.62311667
83 -3.02478333 -7.21067222
84 -5.30039722 -3.02478333
85 -3.36739722 -5.30039722
86 -2.42406389 -3.36739722
87 -2.62450833 -2.42406389
88 -4.42895278 -2.62450833
89 2.43360278 -4.42895278
90 -0.97806389 2.43360278
91 0.41549167 -0.97806389
92 4.09071389 0.41549167
93 3.10015833 4.09071389
94 2.52560278 3.10015833
95 -0.45250833 2.52560278
96 7.12687778 -0.45250833
97 8.55087778 7.12687778
98 9.41821111 8.55087778
99 7.36976667 9.41821111
100 8.37632222 7.36976667
101 2.11187778 8.37632222
102 2.79921111 2.11187778
103 3.02476667 2.79921111
104 2.31098889 3.02476667
105 3.33843333 2.31098889
106 6.06287778 3.33843333
107 5.46376667 6.06287778
108 NA 5.46376667
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -1.17932222 1.73167778
[2,] -1.54798889 -1.17932222
[3,] 1.09256667 -1.54798889
[4,] 5.40712222 1.09256667
[5,] 3.40167778 5.40712222
[6,] 3.17301111 3.40167778
[7,] 1.78556667 3.17301111
[8,] 1.99378889 1.78556667
[9,] -0.66376667 1.99378889
[10,] -1.44532222 -0.66376667
[11,] 1.50856667 -1.44532222
[12,] 2.86295278 1.50856667
[13,] -0.83004722 2.86295278
[14,] 0.57128611 -0.83004722
[15,] 3.04384167 0.57128611
[16,] -1.42060278 3.04384167
[17,] 2.06795278 -1.42060278
[18,] 1.37028611 2.06795278
[19,] -2.48615833 1.37028611
[20,] 0.20506389 -2.48615833
[21,] 0.43450833 0.20506389
[22,] 1.02595278 0.43450833
[23,] 1.97484167 1.02595278
[24,] 3.33122778 1.97484167
[25,] -6.22677222 3.33122778
[26,] -1.30043889 -6.22677222
[27,] 1.93511667 -1.30043889
[28,] -3.33632778 1.93511667
[29,] -2.46677222 -3.33632778
[30,] -0.02043889 -2.46677222
[31,] -2.18988333 -0.02043889
[32,] 2.52633889 -2.18988333
[33,] -4.23821667 2.52633889
[34,] -0.75777222 -4.23821667
[35,] -1.69188333 -0.75777222
[36,] 1.07850278 -1.69188333
[37,] -0.96849722 1.07850278
[38,] 3.73883611 -0.96849722
[39,] -2.66260833 3.73883611
[40,] 0.32094722 -2.66260833
[41,] 0.78750278 0.32094722
[42,] -1.86716389 0.78750278
[43,] -1.73060833 -1.86716389
[44,] -1.70538611 -1.73060833
[45,] -1.42294167 -1.70538611
[46,] 0.81450278 -1.42294167
[47,] 2.12439167 0.81450278
[48,] 1.06777778 2.12439167
[49,] 0.77377778 1.06777778
[50,] 4.17011111 0.77377778
[51,] -0.78633333 4.17011111
[52,] 0.75622222 -0.78633333
[53,] 4.21177778 0.75622222
[54,] 0.09111111 4.21177778
[55,] 5.52166667 0.09111111
[56,] -0.88811111 5.52166667
[57,] 2.27733333 -0.88811111
[58,] 3.69077778 2.27733333
[59,] -1.98233333 3.69077778
[60,] -0.08494722 -1.98233333
[61,] 9.56405278 -0.08494722
[62,] -0.15061389 9.56405278
[63,] 2.21994167 -0.15061389
[64,] 4.62149722 2.21994167
[65,] 0.05605278 4.62149722
[66,] 0.95538611 0.05605278
[67,] 3.74394167 0.95538611
[68,] -0.43083611 3.74394167
[69,] 1.79760833 -0.43083611
[70,] -4.70594722 1.79760833
[71,] -3.92005833 -4.70594722
[72,] -11.81367222 -3.92005833
[73,] -6.31667222 -11.81367222
[74,] -12.47533889 -6.31667222
[75,] -9.58778333 -12.47533889
[76,] -10.29622778 -9.58778333
[77,] -12.60367222 -10.29622778
[78,] -5.52333889 -12.60367222
[79,] -8.08478333 -5.52333889
[80,] -8.10256111 -8.08478333
[81,] -4.62311667 -8.10256111
[82,] -7.21067222 -4.62311667
[83,] -3.02478333 -7.21067222
[84,] -5.30039722 -3.02478333
[85,] -3.36739722 -5.30039722
[86,] -2.42406389 -3.36739722
[87,] -2.62450833 -2.42406389
[88,] -4.42895278 -2.62450833
[89,] 2.43360278 -4.42895278
[90,] -0.97806389 2.43360278
[91,] 0.41549167 -0.97806389
[92,] 4.09071389 0.41549167
[93,] 3.10015833 4.09071389
[94,] 2.52560278 3.10015833
[95,] -0.45250833 2.52560278
[96,] 7.12687778 -0.45250833
[97,] 8.55087778 7.12687778
[98,] 9.41821111 8.55087778
[99,] 7.36976667 9.41821111
[100,] 8.37632222 7.36976667
[101,] 2.11187778 8.37632222
[102,] 2.79921111 2.11187778
[103,] 3.02476667 2.79921111
[104,] 2.31098889 3.02476667
[105,] 3.33843333 2.31098889
[106,] 6.06287778 3.33843333
[107,] 5.46376667 6.06287778
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -1.17932222 1.73167778
2 -1.54798889 -1.17932222
3 1.09256667 -1.54798889
4 5.40712222 1.09256667
5 3.40167778 5.40712222
6 3.17301111 3.40167778
7 1.78556667 3.17301111
8 1.99378889 1.78556667
9 -0.66376667 1.99378889
10 -1.44532222 -0.66376667
11 1.50856667 -1.44532222
12 2.86295278 1.50856667
13 -0.83004722 2.86295278
14 0.57128611 -0.83004722
15 3.04384167 0.57128611
16 -1.42060278 3.04384167
17 2.06795278 -1.42060278
18 1.37028611 2.06795278
19 -2.48615833 1.37028611
20 0.20506389 -2.48615833
21 0.43450833 0.20506389
22 1.02595278 0.43450833
23 1.97484167 1.02595278
24 3.33122778 1.97484167
25 -6.22677222 3.33122778
26 -1.30043889 -6.22677222
27 1.93511667 -1.30043889
28 -3.33632778 1.93511667
29 -2.46677222 -3.33632778
30 -0.02043889 -2.46677222
31 -2.18988333 -0.02043889
32 2.52633889 -2.18988333
33 -4.23821667 2.52633889
34 -0.75777222 -4.23821667
35 -1.69188333 -0.75777222
36 1.07850278 -1.69188333
37 -0.96849722 1.07850278
38 3.73883611 -0.96849722
39 -2.66260833 3.73883611
40 0.32094722 -2.66260833
41 0.78750278 0.32094722
42 -1.86716389 0.78750278
43 -1.73060833 -1.86716389
44 -1.70538611 -1.73060833
45 -1.42294167 -1.70538611
46 0.81450278 -1.42294167
47 2.12439167 0.81450278
48 1.06777778 2.12439167
49 0.77377778 1.06777778
50 4.17011111 0.77377778
51 -0.78633333 4.17011111
52 0.75622222 -0.78633333
53 4.21177778 0.75622222
54 0.09111111 4.21177778
55 5.52166667 0.09111111
56 -0.88811111 5.52166667
57 2.27733333 -0.88811111
58 3.69077778 2.27733333
59 -1.98233333 3.69077778
60 -0.08494722 -1.98233333
61 9.56405278 -0.08494722
62 -0.15061389 9.56405278
63 2.21994167 -0.15061389
64 4.62149722 2.21994167
65 0.05605278 4.62149722
66 0.95538611 0.05605278
67 3.74394167 0.95538611
68 -0.43083611 3.74394167
69 1.79760833 -0.43083611
70 -4.70594722 1.79760833
71 -3.92005833 -4.70594722
72 -11.81367222 -3.92005833
73 -6.31667222 -11.81367222
74 -12.47533889 -6.31667222
75 -9.58778333 -12.47533889
76 -10.29622778 -9.58778333
77 -12.60367222 -10.29622778
78 -5.52333889 -12.60367222
79 -8.08478333 -5.52333889
80 -8.10256111 -8.08478333
81 -4.62311667 -8.10256111
82 -7.21067222 -4.62311667
83 -3.02478333 -7.21067222
84 -5.30039722 -3.02478333
85 -3.36739722 -5.30039722
86 -2.42406389 -3.36739722
87 -2.62450833 -2.42406389
88 -4.42895278 -2.62450833
89 2.43360278 -4.42895278
90 -0.97806389 2.43360278
91 0.41549167 -0.97806389
92 4.09071389 0.41549167
93 3.10015833 4.09071389
94 2.52560278 3.10015833
95 -0.45250833 2.52560278
96 7.12687778 -0.45250833
97 8.55087778 7.12687778
98 9.41821111 8.55087778
99 7.36976667 9.41821111
100 8.37632222 7.36976667
101 2.11187778 8.37632222
102 2.79921111 2.11187778
103 3.02476667 2.79921111
104 2.31098889 3.02476667
105 3.33843333 2.31098889
106 6.06287778 3.33843333
107 5.46376667 6.06287778
> 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/fisher/rcomp/tmp/7ry2q1355321774.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/fisher/rcomp/tmp/8x5191355321774.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/fisher/rcomp/tmp/96ngb1355321774.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/fisher/rcomp/tmp/10jbp01355321774.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/fisher/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/fisher/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/fisher/rcomp/tmp/11uclb1355321774.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/fisher/rcomp/tmp/127og11355321774.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/fisher/rcomp/tmp/13k4871355321774.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/fisher/rcomp/tmp/14atu91355321774.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/fisher/rcomp/tmp/157oqb1355321774.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/fisher/rcomp/tmp/16gno31355321774.tab")
+ }
>
> try(system("convert tmp/1110l1355321774.ps tmp/1110l1355321774.png",intern=TRUE))
character(0)
> try(system("convert tmp/24xq51355321774.ps tmp/24xq51355321774.png",intern=TRUE))
character(0)
> try(system("convert tmp/3c8v61355321774.ps tmp/3c8v61355321774.png",intern=TRUE))
character(0)
> try(system("convert tmp/4s5vj1355321774.ps tmp/4s5vj1355321774.png",intern=TRUE))
character(0)
> try(system("convert tmp/5j15z1355321774.ps tmp/5j15z1355321774.png",intern=TRUE))
character(0)
> try(system("convert tmp/6in0b1355321774.ps tmp/6in0b1355321774.png",intern=TRUE))
character(0)
> try(system("convert tmp/7ry2q1355321774.ps tmp/7ry2q1355321774.png",intern=TRUE))
character(0)
> try(system("convert tmp/8x5191355321774.ps tmp/8x5191355321774.png",intern=TRUE))
character(0)
> try(system("convert tmp/96ngb1355321774.ps tmp/96ngb1355321774.png",intern=TRUE))
character(0)
> try(system("convert tmp/10jbp01355321774.ps tmp/10jbp01355321774.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
6.728 1.555 8.278