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(127,2.75,123,2.75,118,2.55,114,2.5,108,2.5,111,2.1,151,2,159,2,158,2,148,2,138,2,137,2,136,2,133,2,126,2,120,2,114,2,116,2,153,2,162,2,161,2,149,2,139,2,135,2,130,2,127,2,122,2,117,2,112,2,113,2,149,2,157,2,157,2,147,2,137,2,132,2.21,125,2.25,123,2.25,117,2.45,114,2.5,111,2.5,112,2.64,144,2.75,150,2.93,149,3,134,3.17,123,3.25,116,3.39,117,3.5,111,3.5,105,3.65,102,3.75,95,3.75,93,3.9,124,4,130,4,124,4,115,4,106,4,105,4,105,4,101,4,95,4,93,4,84,4,87,4,116,4.18,120,4.25,117,4.25,109,3.97,105,3.42,107,2.75),dim=c(2,72),dimnames=list(c('Werkloosheid','Rente'),1:72))
> y <- array(NA,dim=c(2,72),dimnames=list(c('Werkloosheid','Rente'),1:72))
> 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
Werkloosheid Rente M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 127 2.75 1 0 0 0 0 0 0 0 0 0 0 1
2 123 2.75 0 1 0 0 0 0 0 0 0 0 0 2
3 118 2.55 0 0 1 0 0 0 0 0 0 0 0 3
4 114 2.50 0 0 0 1 0 0 0 0 0 0 0 4
5 108 2.50 0 0 0 0 1 0 0 0 0 0 0 5
6 111 2.10 0 0 0 0 0 1 0 0 0 0 0 6
7 151 2.00 0 0 0 0 0 0 1 0 0 0 0 7
8 159 2.00 0 0 0 0 0 0 0 1 0 0 0 8
9 158 2.00 0 0 0 0 0 0 0 0 1 0 0 9
10 148 2.00 0 0 0 0 0 0 0 0 0 1 0 10
11 138 2.00 0 0 0 0 0 0 0 0 0 0 1 11
12 137 2.00 0 0 0 0 0 0 0 0 0 0 0 12
13 136 2.00 1 0 0 0 0 0 0 0 0 0 0 13
14 133 2.00 0 1 0 0 0 0 0 0 0 0 0 14
15 126 2.00 0 0 1 0 0 0 0 0 0 0 0 15
16 120 2.00 0 0 0 1 0 0 0 0 0 0 0 16
17 114 2.00 0 0 0 0 1 0 0 0 0 0 0 17
18 116 2.00 0 0 0 0 0 1 0 0 0 0 0 18
19 153 2.00 0 0 0 0 0 0 1 0 0 0 0 19
20 162 2.00 0 0 0 0 0 0 0 1 0 0 0 20
21 161 2.00 0 0 0 0 0 0 0 0 1 0 0 21
22 149 2.00 0 0 0 0 0 0 0 0 0 1 0 22
23 139 2.00 0 0 0 0 0 0 0 0 0 0 1 23
24 135 2.00 0 0 0 0 0 0 0 0 0 0 0 24
25 130 2.00 1 0 0 0 0 0 0 0 0 0 0 25
26 127 2.00 0 1 0 0 0 0 0 0 0 0 0 26
27 122 2.00 0 0 1 0 0 0 0 0 0 0 0 27
28 117 2.00 0 0 0 1 0 0 0 0 0 0 0 28
29 112 2.00 0 0 0 0 1 0 0 0 0 0 0 29
30 113 2.00 0 0 0 0 0 1 0 0 0 0 0 30
31 149 2.00 0 0 0 0 0 0 1 0 0 0 0 31
32 157 2.00 0 0 0 0 0 0 0 1 0 0 0 32
33 157 2.00 0 0 0 0 0 0 0 0 1 0 0 33
34 147 2.00 0 0 0 0 0 0 0 0 0 1 0 34
35 137 2.00 0 0 0 0 0 0 0 0 0 0 1 35
36 132 2.21 0 0 0 0 0 0 0 0 0 0 0 36
37 125 2.25 1 0 0 0 0 0 0 0 0 0 0 37
38 123 2.25 0 1 0 0 0 0 0 0 0 0 0 38
39 117 2.45 0 0 1 0 0 0 0 0 0 0 0 39
40 114 2.50 0 0 0 1 0 0 0 0 0 0 0 40
41 111 2.50 0 0 0 0 1 0 0 0 0 0 0 41
42 112 2.64 0 0 0 0 0 1 0 0 0 0 0 42
43 144 2.75 0 0 0 0 0 0 1 0 0 0 0 43
44 150 2.93 0 0 0 0 0 0 0 1 0 0 0 44
45 149 3.00 0 0 0 0 0 0 0 0 1 0 0 45
46 134 3.17 0 0 0 0 0 0 0 0 0 1 0 46
47 123 3.25 0 0 0 0 0 0 0 0 0 0 1 47
48 116 3.39 0 0 0 0 0 0 0 0 0 0 0 48
49 117 3.50 1 0 0 0 0 0 0 0 0 0 0 49
50 111 3.50 0 1 0 0 0 0 0 0 0 0 0 50
51 105 3.65 0 0 1 0 0 0 0 0 0 0 0 51
52 102 3.75 0 0 0 1 0 0 0 0 0 0 0 52
53 95 3.75 0 0 0 0 1 0 0 0 0 0 0 53
54 93 3.90 0 0 0 0 0 1 0 0 0 0 0 54
55 124 4.00 0 0 0 0 0 0 1 0 0 0 0 55
56 130 4.00 0 0 0 0 0 0 0 1 0 0 0 56
57 124 4.00 0 0 0 0 0 0 0 0 1 0 0 57
58 115 4.00 0 0 0 0 0 0 0 0 0 1 0 58
59 106 4.00 0 0 0 0 0 0 0 0 0 0 1 59
60 105 4.00 0 0 0 0 0 0 0 0 0 0 0 60
61 105 4.00 1 0 0 0 0 0 0 0 0 0 0 61
62 101 4.00 0 1 0 0 0 0 0 0 0 0 0 62
63 95 4.00 0 0 1 0 0 0 0 0 0 0 0 63
64 93 4.00 0 0 0 1 0 0 0 0 0 0 0 64
65 84 4.00 0 0 0 0 1 0 0 0 0 0 0 65
66 87 4.00 0 0 0 0 0 1 0 0 0 0 0 66
67 116 4.18 0 0 0 0 0 0 1 0 0 0 0 67
68 120 4.25 0 0 0 0 0 0 0 1 0 0 0 68
69 117 4.25 0 0 0 0 0 0 0 0 1 0 0 69
70 109 3.97 0 0 0 0 0 0 0 0 0 1 0 70
71 105 3.42 0 0 0 0 0 0 0 0 0 0 1 71
72 107 2.75 0 0 0 0 0 0 0 0 0 0 0 72
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Rente M1 M2 M3 M4
160.40921 -11.78875 -0.01799 -3.53501 -8.92399 -12.41120
M5 M6 M7 M8 M9 M10
-18.26156 -16.99472 17.89138 25.36555 23.65273 12.91957
M11 t
3.14576 -0.14964
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-10.2161 -2.2599 0.1651 2.4834 7.0381
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 160.40921 2.27414 70.536 < 2e-16 ***
Rente -11.78875 1.00994 -11.673 < 2e-16 ***
M1 -0.01799 2.32665 -0.008 0.993859
M2 -3.53501 2.31849 -1.525 0.132768
M3 -8.92399 2.31488 -3.855 0.000292 ***
M4 -12.41120 2.31024 -5.372 1.44e-06 ***
M5 -18.26156 2.30363 -7.927 8.11e-11 ***
M6 -16.99472 2.29561 -7.403 6.16e-10 ***
M7 17.89138 2.29599 7.792 1.37e-10 ***
M8 25.36555 2.29586 11.048 6.78e-16 ***
M9 23.65273 2.29241 10.318 9.61e-15 ***
M10 12.91957 2.28652 5.650 5.10e-07 ***
M11 3.14576 2.27892 1.380 0.172769
t -0.14964 0.04149 -3.607 0.000646 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3.944 on 58 degrees of freedom
Multiple R-squared: 0.9656, Adjusted R-squared: 0.9579
F-statistic: 125.3 on 13 and 58 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,] 4.580512e-03 9.161024e-03 0.995419488
[2,] 4.053346e-03 8.106691e-03 0.995946654
[3,] 1.274063e-03 2.548126e-03 0.998725937
[4,] 2.377394e-04 4.754788e-04 0.999762261
[5,] 4.046086e-05 8.092171e-05 0.999959539
[6,] 2.921423e-05 5.842846e-05 0.999970786
[7,] 1.225136e-05 2.450272e-05 0.999987749
[8,] 1.691598e-04 3.383197e-04 0.999830840
[9,] 1.953594e-02 3.907188e-02 0.980464058
[10,] 4.698965e-02 9.397929e-02 0.953010355
[11,] 4.349731e-02 8.699462e-02 0.956502688
[12,] 5.496178e-02 1.099236e-01 0.945038219
[13,] 6.323198e-02 1.264640e-01 0.936768024
[14,] 9.820163e-02 1.964033e-01 0.901798367
[15,] 9.916514e-02 1.983303e-01 0.900834858
[16,] 8.690241e-02 1.738048e-01 0.913097586
[17,] 5.754428e-02 1.150886e-01 0.942455718
[18,] 3.631018e-02 7.262036e-02 0.963689818
[19,] 2.190609e-02 4.381218e-02 0.978093908
[20,] 1.281558e-02 2.563115e-02 0.987184424
[21,] 5.254505e-02 1.050901e-01 0.947454953
[22,] 8.850714e-02 1.770143e-01 0.911492857
[23,] 1.528856e-01 3.057713e-01 0.847114350
[24,] 5.420209e-01 9.159582e-01 0.457979114
[25,] 7.471710e-01 5.056580e-01 0.252829000
[26,] 8.701945e-01 2.596111e-01 0.129805547
[27,] 8.450436e-01 3.099128e-01 0.154956415
[28,] 7.825310e-01 4.349379e-01 0.217468962
[29,] 9.565833e-01 8.683340e-02 0.043416698
[30,] 9.891904e-01 2.161911e-02 0.010809557
[31,] 9.983140e-01 3.372009e-03 0.001686004
[32,] 9.970832e-01 5.833551e-03 0.002916775
[33,] 9.952346e-01 9.530826e-03 0.004765413
[34,] 9.893151e-01 2.136971e-02 0.010684854
[35,] 9.757816e-01 4.843686e-02 0.024218430
[36,] 9.496678e-01 1.006644e-01 0.050332218
[37,] 9.721982e-01 5.560354e-02 0.027801771
[38,] 9.450625e-01 1.098750e-01 0.054937503
[39,] 8.807097e-01 2.385807e-01 0.119290346
> postscript(file="/var/www/html/rcomp/tmp/1o0481258711308.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/22asl1258711308.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/34p7n1258711308.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/4agnq1258711308.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/5tswy1258711308.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 = 72
Frequency = 1
1 2 3 4 5 6
-0.82252793 -1.15586127 -2.97499632 -3.92757954 -3.92757954 -6.76028518
7 8 9 10 11 12
-2.67561619 -2.00014737 -1.13768277 -0.25488905 -0.33143709 1.96396282
13 14 15 16 17 18
1.13158862 1.79825529 0.33686991 -2.02627590 -2.02627590 -1.14348218
19 20 21 22 23 24
1.12006165 2.79553047 3.65799507 2.54078879 2.46424075 1.75964066
25 26 27 28 29 30
-3.07273354 -2.40606687 -1.86745225 -3.23059806 -2.23059806 -2.34780434
31 32 33 34 35 36
-1.08426051 -0.40879169 1.45367291 2.33646663 2.25991859 3.03095566
37 38 39 40 41 42
-3.32986861 -1.66320194 0.23316236 1.45945397 4.45945397 5.99267247
43 44 45 46 47 48
4.55297862 5.35042215 7.03809913 4.92498008 4.79153191 2.73735660
49 50 51 52 53 54
5.20174471 2.86841138 4.17533826 5.99106729 4.99106729 3.64217327
55 56 57 58 59 60
1.08459194 -0.23993924 -4.37747464 -2.49468092 -1.57122896 0.72417095
61 62 63 64 65 66
0.89179675 0.55846341 0.09707804 1.73393223 -1.26606777 0.61672595
67 68 69 70 71 72
-2.99775551 -5.49707431 -6.63460971 -7.05266553 -7.61302519 -10.21608669
> postscript(file="/var/www/html/rcomp/tmp/6uue31258711308.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 = 72
Frequency = 1
lag(myerror, k = 1) myerror
0 -0.82252793 NA
1 -1.15586127 -0.82252793
2 -2.97499632 -1.15586127
3 -3.92757954 -2.97499632
4 -3.92757954 -3.92757954
5 -6.76028518 -3.92757954
6 -2.67561619 -6.76028518
7 -2.00014737 -2.67561619
8 -1.13768277 -2.00014737
9 -0.25488905 -1.13768277
10 -0.33143709 -0.25488905
11 1.96396282 -0.33143709
12 1.13158862 1.96396282
13 1.79825529 1.13158862
14 0.33686991 1.79825529
15 -2.02627590 0.33686991
16 -2.02627590 -2.02627590
17 -1.14348218 -2.02627590
18 1.12006165 -1.14348218
19 2.79553047 1.12006165
20 3.65799507 2.79553047
21 2.54078879 3.65799507
22 2.46424075 2.54078879
23 1.75964066 2.46424075
24 -3.07273354 1.75964066
25 -2.40606687 -3.07273354
26 -1.86745225 -2.40606687
27 -3.23059806 -1.86745225
28 -2.23059806 -3.23059806
29 -2.34780434 -2.23059806
30 -1.08426051 -2.34780434
31 -0.40879169 -1.08426051
32 1.45367291 -0.40879169
33 2.33646663 1.45367291
34 2.25991859 2.33646663
35 3.03095566 2.25991859
36 -3.32986861 3.03095566
37 -1.66320194 -3.32986861
38 0.23316236 -1.66320194
39 1.45945397 0.23316236
40 4.45945397 1.45945397
41 5.99267247 4.45945397
42 4.55297862 5.99267247
43 5.35042215 4.55297862
44 7.03809913 5.35042215
45 4.92498008 7.03809913
46 4.79153191 4.92498008
47 2.73735660 4.79153191
48 5.20174471 2.73735660
49 2.86841138 5.20174471
50 4.17533826 2.86841138
51 5.99106729 4.17533826
52 4.99106729 5.99106729
53 3.64217327 4.99106729
54 1.08459194 3.64217327
55 -0.23993924 1.08459194
56 -4.37747464 -0.23993924
57 -2.49468092 -4.37747464
58 -1.57122896 -2.49468092
59 0.72417095 -1.57122896
60 0.89179675 0.72417095
61 0.55846341 0.89179675
62 0.09707804 0.55846341
63 1.73393223 0.09707804
64 -1.26606777 1.73393223
65 0.61672595 -1.26606777
66 -2.99775551 0.61672595
67 -5.49707431 -2.99775551
68 -6.63460971 -5.49707431
69 -7.05266553 -6.63460971
70 -7.61302519 -7.05266553
71 -10.21608669 -7.61302519
72 NA -10.21608669
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -1.15586127 -0.82252793
[2,] -2.97499632 -1.15586127
[3,] -3.92757954 -2.97499632
[4,] -3.92757954 -3.92757954
[5,] -6.76028518 -3.92757954
[6,] -2.67561619 -6.76028518
[7,] -2.00014737 -2.67561619
[8,] -1.13768277 -2.00014737
[9,] -0.25488905 -1.13768277
[10,] -0.33143709 -0.25488905
[11,] 1.96396282 -0.33143709
[12,] 1.13158862 1.96396282
[13,] 1.79825529 1.13158862
[14,] 0.33686991 1.79825529
[15,] -2.02627590 0.33686991
[16,] -2.02627590 -2.02627590
[17,] -1.14348218 -2.02627590
[18,] 1.12006165 -1.14348218
[19,] 2.79553047 1.12006165
[20,] 3.65799507 2.79553047
[21,] 2.54078879 3.65799507
[22,] 2.46424075 2.54078879
[23,] 1.75964066 2.46424075
[24,] -3.07273354 1.75964066
[25,] -2.40606687 -3.07273354
[26,] -1.86745225 -2.40606687
[27,] -3.23059806 -1.86745225
[28,] -2.23059806 -3.23059806
[29,] -2.34780434 -2.23059806
[30,] -1.08426051 -2.34780434
[31,] -0.40879169 -1.08426051
[32,] 1.45367291 -0.40879169
[33,] 2.33646663 1.45367291
[34,] 2.25991859 2.33646663
[35,] 3.03095566 2.25991859
[36,] -3.32986861 3.03095566
[37,] -1.66320194 -3.32986861
[38,] 0.23316236 -1.66320194
[39,] 1.45945397 0.23316236
[40,] 4.45945397 1.45945397
[41,] 5.99267247 4.45945397
[42,] 4.55297862 5.99267247
[43,] 5.35042215 4.55297862
[44,] 7.03809913 5.35042215
[45,] 4.92498008 7.03809913
[46,] 4.79153191 4.92498008
[47,] 2.73735660 4.79153191
[48,] 5.20174471 2.73735660
[49,] 2.86841138 5.20174471
[50,] 4.17533826 2.86841138
[51,] 5.99106729 4.17533826
[52,] 4.99106729 5.99106729
[53,] 3.64217327 4.99106729
[54,] 1.08459194 3.64217327
[55,] -0.23993924 1.08459194
[56,] -4.37747464 -0.23993924
[57,] -2.49468092 -4.37747464
[58,] -1.57122896 -2.49468092
[59,] 0.72417095 -1.57122896
[60,] 0.89179675 0.72417095
[61,] 0.55846341 0.89179675
[62,] 0.09707804 0.55846341
[63,] 1.73393223 0.09707804
[64,] -1.26606777 1.73393223
[65,] 0.61672595 -1.26606777
[66,] -2.99775551 0.61672595
[67,] -5.49707431 -2.99775551
[68,] -6.63460971 -5.49707431
[69,] -7.05266553 -6.63460971
[70,] -7.61302519 -7.05266553
[71,] -10.21608669 -7.61302519
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -1.15586127 -0.82252793
2 -2.97499632 -1.15586127
3 -3.92757954 -2.97499632
4 -3.92757954 -3.92757954
5 -6.76028518 -3.92757954
6 -2.67561619 -6.76028518
7 -2.00014737 -2.67561619
8 -1.13768277 -2.00014737
9 -0.25488905 -1.13768277
10 -0.33143709 -0.25488905
11 1.96396282 -0.33143709
12 1.13158862 1.96396282
13 1.79825529 1.13158862
14 0.33686991 1.79825529
15 -2.02627590 0.33686991
16 -2.02627590 -2.02627590
17 -1.14348218 -2.02627590
18 1.12006165 -1.14348218
19 2.79553047 1.12006165
20 3.65799507 2.79553047
21 2.54078879 3.65799507
22 2.46424075 2.54078879
23 1.75964066 2.46424075
24 -3.07273354 1.75964066
25 -2.40606687 -3.07273354
26 -1.86745225 -2.40606687
27 -3.23059806 -1.86745225
28 -2.23059806 -3.23059806
29 -2.34780434 -2.23059806
30 -1.08426051 -2.34780434
31 -0.40879169 -1.08426051
32 1.45367291 -0.40879169
33 2.33646663 1.45367291
34 2.25991859 2.33646663
35 3.03095566 2.25991859
36 -3.32986861 3.03095566
37 -1.66320194 -3.32986861
38 0.23316236 -1.66320194
39 1.45945397 0.23316236
40 4.45945397 1.45945397
41 5.99267247 4.45945397
42 4.55297862 5.99267247
43 5.35042215 4.55297862
44 7.03809913 5.35042215
45 4.92498008 7.03809913
46 4.79153191 4.92498008
47 2.73735660 4.79153191
48 5.20174471 2.73735660
49 2.86841138 5.20174471
50 4.17533826 2.86841138
51 5.99106729 4.17533826
52 4.99106729 5.99106729
53 3.64217327 4.99106729
54 1.08459194 3.64217327
55 -0.23993924 1.08459194
56 -4.37747464 -0.23993924
57 -2.49468092 -4.37747464
58 -1.57122896 -2.49468092
59 0.72417095 -1.57122896
60 0.89179675 0.72417095
61 0.55846341 0.89179675
62 0.09707804 0.55846341
63 1.73393223 0.09707804
64 -1.26606777 1.73393223
65 0.61672595 -1.26606777
66 -2.99775551 0.61672595
67 -5.49707431 -2.99775551
68 -6.63460971 -5.49707431
69 -7.05266553 -6.63460971
70 -7.61302519 -7.05266553
71 -10.21608669 -7.61302519
> 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/7bd321258711308.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/8oh2u1258711308.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/9o2pu1258711308.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/10atl61258711308.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/11slot1258711308.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/12kezz1258711308.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/13vm561258711308.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/14tfxb1258711308.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/15fik81258711308.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/16fmat1258711308.tab")
+ }
>
> system("convert tmp/1o0481258711308.ps tmp/1o0481258711308.png")
> system("convert tmp/22asl1258711308.ps tmp/22asl1258711308.png")
> system("convert tmp/34p7n1258711308.ps tmp/34p7n1258711308.png")
> system("convert tmp/4agnq1258711308.ps tmp/4agnq1258711308.png")
> system("convert tmp/5tswy1258711308.ps tmp/5tswy1258711308.png")
> system("convert tmp/6uue31258711308.ps tmp/6uue31258711308.png")
> system("convert tmp/7bd321258711308.ps tmp/7bd321258711308.png")
> system("convert tmp/8oh2u1258711308.ps tmp/8oh2u1258711308.png")
> system("convert tmp/9o2pu1258711308.ps tmp/9o2pu1258711308.png")
> system("convert tmp/10atl61258711308.ps tmp/10atl61258711308.png")
>
>
> proc.time()
user system elapsed
2.512 1.543 3.252