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 = 'No 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
1 127 2.75 1 0 0 0 0 0 0 0 0 0 0
2 123 2.75 0 1 0 0 0 0 0 0 0 0 0
3 118 2.55 0 0 1 0 0 0 0 0 0 0 0
4 114 2.50 0 0 0 1 0 0 0 0 0 0 0
5 108 2.50 0 0 0 0 1 0 0 0 0 0 0
6 111 2.10 0 0 0 0 0 1 0 0 0 0 0
7 151 2.00 0 0 0 0 0 0 1 0 0 0 0
8 159 2.00 0 0 0 0 0 0 0 1 0 0 0
9 158 2.00 0 0 0 0 0 0 0 0 1 0 0
10 148 2.00 0 0 0 0 0 0 0 0 0 1 0
11 138 2.00 0 0 0 0 0 0 0 0 0 0 1
12 137 2.00 0 0 0 0 0 0 0 0 0 0 0
13 136 2.00 1 0 0 0 0 0 0 0 0 0 0
14 133 2.00 0 1 0 0 0 0 0 0 0 0 0
15 126 2.00 0 0 1 0 0 0 0 0 0 0 0
16 120 2.00 0 0 0 1 0 0 0 0 0 0 0
17 114 2.00 0 0 0 0 1 0 0 0 0 0 0
18 116 2.00 0 0 0 0 0 1 0 0 0 0 0
19 153 2.00 0 0 0 0 0 0 1 0 0 0 0
20 162 2.00 0 0 0 0 0 0 0 1 0 0 0
21 161 2.00 0 0 0 0 0 0 0 0 1 0 0
22 149 2.00 0 0 0 0 0 0 0 0 0 1 0
23 139 2.00 0 0 0 0 0 0 0 0 0 0 1
24 135 2.00 0 0 0 0 0 0 0 0 0 0 0
25 130 2.00 1 0 0 0 0 0 0 0 0 0 0
26 127 2.00 0 1 0 0 0 0 0 0 0 0 0
27 122 2.00 0 0 1 0 0 0 0 0 0 0 0
28 117 2.00 0 0 0 1 0 0 0 0 0 0 0
29 112 2.00 0 0 0 0 1 0 0 0 0 0 0
30 113 2.00 0 0 0 0 0 1 0 0 0 0 0
31 149 2.00 0 0 0 0 0 0 1 0 0 0 0
32 157 2.00 0 0 0 0 0 0 0 1 0 0 0
33 157 2.00 0 0 0 0 0 0 0 0 1 0 0
34 147 2.00 0 0 0 0 0 0 0 0 0 1 0
35 137 2.00 0 0 0 0 0 0 0 0 0 0 1
36 132 2.21 0 0 0 0 0 0 0 0 0 0 0
37 125 2.25 1 0 0 0 0 0 0 0 0 0 0
38 123 2.25 0 1 0 0 0 0 0 0 0 0 0
39 117 2.45 0 0 1 0 0 0 0 0 0 0 0
40 114 2.50 0 0 0 1 0 0 0 0 0 0 0
41 111 2.50 0 0 0 0 1 0 0 0 0 0 0
42 112 2.64 0 0 0 0 0 1 0 0 0 0 0
43 144 2.75 0 0 0 0 0 0 1 0 0 0 0
44 150 2.93 0 0 0 0 0 0 0 1 0 0 0
45 149 3.00 0 0 0 0 0 0 0 0 1 0 0
46 134 3.17 0 0 0 0 0 0 0 0 0 1 0
47 123 3.25 0 0 0 0 0 0 0 0 0 0 1
48 116 3.39 0 0 0 0 0 0 0 0 0 0 0
49 117 3.50 1 0 0 0 0 0 0 0 0 0 0
50 111 3.50 0 1 0 0 0 0 0 0 0 0 0
51 105 3.65 0 0 1 0 0 0 0 0 0 0 0
52 102 3.75 0 0 0 1 0 0 0 0 0 0 0
53 95 3.75 0 0 0 0 1 0 0 0 0 0 0
54 93 3.90 0 0 0 0 0 1 0 0 0 0 0
55 124 4.00 0 0 0 0 0 0 1 0 0 0 0
56 130 4.00 0 0 0 0 0 0 0 1 0 0 0
57 124 4.00 0 0 0 0 0 0 0 0 1 0 0
58 115 4.00 0 0 0 0 0 0 0 0 0 1 0
59 106 4.00 0 0 0 0 0 0 0 0 0 0 1
60 105 4.00 0 0 0 0 0 0 0 0 0 0 0
61 105 4.00 1 0 0 0 0 0 0 0 0 0 0
62 101 4.00 0 1 0 0 0 0 0 0 0 0 0
63 95 4.00 0 0 1 0 0 0 0 0 0 0 0
64 93 4.00 0 0 0 1 0 0 0 0 0 0 0
65 84 4.00 0 0 0 0 1 0 0 0 0 0 0
66 87 4.00 0 0 0 0 0 1 0 0 0 0 0
67 116 4.18 0 0 0 0 0 0 1 0 0 0 0
68 120 4.25 0 0 0 0 0 0 0 1 0 0 0
69 117 4.25 0 0 0 0 0 0 0 0 1 0 0
70 109 3.97 0 0 0 0 0 0 0 0 0 1 0
71 105 3.42 0 0 0 0 0 0 0 0 0 0 1
72 107 2.75 0 0 0 0 0 0 0 0 0 0 0
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Rente M1 M2 M3 M4
162.436 -14.839 1.704 -1.962 -7.425 -11.011
M5 M6 M7 M8 M9 M10
-17.011 -15.949 18.934 26.386 24.559 13.620
M11
3.458
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-14.6290 -2.0883 0.5774 2.6999 6.5215
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 162.4360 2.4175 67.192 < 2e-16 ***
Rente -14.8389 0.6057 -24.500 < 2e-16 ***
M1 1.7043 2.4981 0.682 0.49776
M2 -1.9624 2.4981 -0.786 0.43529
M3 -7.4247 2.4983 -2.972 0.00428 **
M4 -11.0107 2.4984 -4.407 4.49e-05 ***
M5 -17.0107 2.4984 -6.809 5.70e-09 ***
M6 -15.9495 2.4983 -6.384 2.95e-08 ***
M7 18.9344 2.4988 7.577 2.84e-10 ***
M8 26.3860 2.4995 10.557 3.25e-15 ***
M9 24.5592 2.4997 9.825 4.97e-14 ***
M10 13.6205 2.4994 5.450 1.04e-06 ***
M11 3.4581 2.4983 1.384 0.17152
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 4.327 on 59 degrees of freedom
Multiple R-squared: 0.9579, Adjusted R-squared: 0.9494
F-statistic: 111.9 on 12 and 59 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,] 1.682698e-03 3.365397e-03 0.9983173
[2,] 1.784031e-04 3.568063e-04 0.9998216
[3,] 2.986826e-03 5.973652e-03 0.9970132
[4,] 1.249093e-03 2.498185e-03 0.9987509
[5,] 1.007247e-03 2.014493e-03 0.9989928
[6,] 7.568044e-04 1.513609e-03 0.9992432
[7,] 2.442006e-04 4.884012e-04 0.9997558
[8,] 7.781770e-05 1.556354e-04 0.9999222
[9,] 3.500213e-05 7.000426e-05 0.9999650
[10,] 5.213604e-04 1.042721e-03 0.9994786
[11,] 1.018963e-03 2.037927e-03 0.9989810
[12,] 6.428701e-04 1.285740e-03 0.9993571
[13,] 4.314030e-04 8.628060e-04 0.9995686
[14,] 2.201048e-04 4.402095e-04 0.9997799
[15,] 1.193096e-04 2.386192e-04 0.9998807
[16,] 8.582385e-05 1.716477e-04 0.9999142
[17,] 6.976699e-05 1.395340e-04 0.9999302
[18,] 3.723181e-05 7.446362e-05 0.9999628
[19,] 1.564308e-05 3.128617e-05 0.9999844
[20,] 6.473295e-06 1.294659e-05 0.9999935
[21,] 3.513073e-06 7.026147e-06 0.9999965
[22,] 2.198725e-05 4.397450e-05 0.9999780
[23,] 3.936544e-05 7.873087e-05 0.9999606
[24,] 2.853344e-05 5.706689e-05 0.9999715
[25,] 1.976562e-05 3.953124e-05 0.9999802
[26,] 1.289084e-05 2.578168e-05 0.9999871
[27,] 1.104073e-05 2.208147e-05 0.9999890
[28,] 4.282431e-06 8.564862e-06 0.9999957
[29,] 1.772124e-06 3.544247e-06 0.9999982
[30,] 1.896380e-06 3.792761e-06 0.9999981
[31,] 5.281333e-06 1.056267e-05 0.9999947
[32,] 3.614741e-05 7.229483e-05 0.9999639
[33,] 2.547260e-04 5.094521e-04 0.9997453
[34,] 2.730806e-04 5.461612e-04 0.9997269
[35,] 2.483952e-04 4.967904e-04 0.9997516
[36,] 3.165298e-04 6.330596e-04 0.9996835
[37,] 5.014186e-04 1.002837e-03 0.9994986
[38,] 2.137219e-03 4.274439e-03 0.9978628
[39,] 2.171991e-03 4.343982e-03 0.9978280
[40,] 7.643852e-03 1.528770e-02 0.9923561
[41,] 9.018269e-02 1.803654e-01 0.9098173
> postscript(file="/var/www/html/rcomp/tmp/1bysl1258711054.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/2xj5c1258711054.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/3900q1258711054.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/4a9jj1258711054.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/5h2s81258711054.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
3.6666667 3.3333333 0.8279142 -0.3280125 -0.3280125 -4.3248593
7 8 9 10 11 12
-0.6926295 -0.1442503 0.6826292 1.6213424 1.7837229 4.2417975
13 14 15 16 17 18
1.5374917 2.2041584 0.6665192 -1.7474625 -1.7474625 -0.8087493
19 20 21 22 23 24
1.3073705 2.8557497 3.6826292 2.6213424 2.7837229 2.2417975
25 26 27 28 29 30
-4.4625083 -3.7958416 -3.3334808 -4.7474625 -3.7474625 -3.8087493
31 32 33 34 35 36
-2.6926295 -2.1442503 -0.3173708 0.6213424 0.7837229 2.3579665
37 38 39 40 41 42
-5.7527833 -4.0861166 -1.6559758 -0.3280125 2.6719875 4.6881467
43 44 45 46 47 48
3.4365455 4.6559267 6.5215292 4.9828553 5.3323478 3.8678685
49 50 51 52 53 54
4.7958416 2.4625083 4.1507041 6.2206124 5.2206124 4.3851606
55 56 57 58 59 60
1.9851704 0.5335496 -3.6395709 -1.7008577 -0.5384772 1.9195974
61 62 63 64 65 66
0.2152916 -0.1180417 -0.6556809 0.9303374 -2.0696626 -0.1309494
67 68 69 70 71 72
-3.3438276 -5.7567254 -6.9298459 -8.1460247 -10.1450392 -14.6290275
> postscript(file="/var/www/html/rcomp/tmp/627bg1258711054.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 3.6666667 NA
1 3.3333333 3.6666667
2 0.8279142 3.3333333
3 -0.3280125 0.8279142
4 -0.3280125 -0.3280125
5 -4.3248593 -0.3280125
6 -0.6926295 -4.3248593
7 -0.1442503 -0.6926295
8 0.6826292 -0.1442503
9 1.6213424 0.6826292
10 1.7837229 1.6213424
11 4.2417975 1.7837229
12 1.5374917 4.2417975
13 2.2041584 1.5374917
14 0.6665192 2.2041584
15 -1.7474625 0.6665192
16 -1.7474625 -1.7474625
17 -0.8087493 -1.7474625
18 1.3073705 -0.8087493
19 2.8557497 1.3073705
20 3.6826292 2.8557497
21 2.6213424 3.6826292
22 2.7837229 2.6213424
23 2.2417975 2.7837229
24 -4.4625083 2.2417975
25 -3.7958416 -4.4625083
26 -3.3334808 -3.7958416
27 -4.7474625 -3.3334808
28 -3.7474625 -4.7474625
29 -3.8087493 -3.7474625
30 -2.6926295 -3.8087493
31 -2.1442503 -2.6926295
32 -0.3173708 -2.1442503
33 0.6213424 -0.3173708
34 0.7837229 0.6213424
35 2.3579665 0.7837229
36 -5.7527833 2.3579665
37 -4.0861166 -5.7527833
38 -1.6559758 -4.0861166
39 -0.3280125 -1.6559758
40 2.6719875 -0.3280125
41 4.6881467 2.6719875
42 3.4365455 4.6881467
43 4.6559267 3.4365455
44 6.5215292 4.6559267
45 4.9828553 6.5215292
46 5.3323478 4.9828553
47 3.8678685 5.3323478
48 4.7958416 3.8678685
49 2.4625083 4.7958416
50 4.1507041 2.4625083
51 6.2206124 4.1507041
52 5.2206124 6.2206124
53 4.3851606 5.2206124
54 1.9851704 4.3851606
55 0.5335496 1.9851704
56 -3.6395709 0.5335496
57 -1.7008577 -3.6395709
58 -0.5384772 -1.7008577
59 1.9195974 -0.5384772
60 0.2152916 1.9195974
61 -0.1180417 0.2152916
62 -0.6556809 -0.1180417
63 0.9303374 -0.6556809
64 -2.0696626 0.9303374
65 -0.1309494 -2.0696626
66 -3.3438276 -0.1309494
67 -5.7567254 -3.3438276
68 -6.9298459 -5.7567254
69 -8.1460247 -6.9298459
70 -10.1450392 -8.1460247
71 -14.6290275 -10.1450392
72 NA -14.6290275
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 3.3333333 3.6666667
[2,] 0.8279142 3.3333333
[3,] -0.3280125 0.8279142
[4,] -0.3280125 -0.3280125
[5,] -4.3248593 -0.3280125
[6,] -0.6926295 -4.3248593
[7,] -0.1442503 -0.6926295
[8,] 0.6826292 -0.1442503
[9,] 1.6213424 0.6826292
[10,] 1.7837229 1.6213424
[11,] 4.2417975 1.7837229
[12,] 1.5374917 4.2417975
[13,] 2.2041584 1.5374917
[14,] 0.6665192 2.2041584
[15,] -1.7474625 0.6665192
[16,] -1.7474625 -1.7474625
[17,] -0.8087493 -1.7474625
[18,] 1.3073705 -0.8087493
[19,] 2.8557497 1.3073705
[20,] 3.6826292 2.8557497
[21,] 2.6213424 3.6826292
[22,] 2.7837229 2.6213424
[23,] 2.2417975 2.7837229
[24,] -4.4625083 2.2417975
[25,] -3.7958416 -4.4625083
[26,] -3.3334808 -3.7958416
[27,] -4.7474625 -3.3334808
[28,] -3.7474625 -4.7474625
[29,] -3.8087493 -3.7474625
[30,] -2.6926295 -3.8087493
[31,] -2.1442503 -2.6926295
[32,] -0.3173708 -2.1442503
[33,] 0.6213424 -0.3173708
[34,] 0.7837229 0.6213424
[35,] 2.3579665 0.7837229
[36,] -5.7527833 2.3579665
[37,] -4.0861166 -5.7527833
[38,] -1.6559758 -4.0861166
[39,] -0.3280125 -1.6559758
[40,] 2.6719875 -0.3280125
[41,] 4.6881467 2.6719875
[42,] 3.4365455 4.6881467
[43,] 4.6559267 3.4365455
[44,] 6.5215292 4.6559267
[45,] 4.9828553 6.5215292
[46,] 5.3323478 4.9828553
[47,] 3.8678685 5.3323478
[48,] 4.7958416 3.8678685
[49,] 2.4625083 4.7958416
[50,] 4.1507041 2.4625083
[51,] 6.2206124 4.1507041
[52,] 5.2206124 6.2206124
[53,] 4.3851606 5.2206124
[54,] 1.9851704 4.3851606
[55,] 0.5335496 1.9851704
[56,] -3.6395709 0.5335496
[57,] -1.7008577 -3.6395709
[58,] -0.5384772 -1.7008577
[59,] 1.9195974 -0.5384772
[60,] 0.2152916 1.9195974
[61,] -0.1180417 0.2152916
[62,] -0.6556809 -0.1180417
[63,] 0.9303374 -0.6556809
[64,] -2.0696626 0.9303374
[65,] -0.1309494 -2.0696626
[66,] -3.3438276 -0.1309494
[67,] -5.7567254 -3.3438276
[68,] -6.9298459 -5.7567254
[69,] -8.1460247 -6.9298459
[70,] -10.1450392 -8.1460247
[71,] -14.6290275 -10.1450392
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 3.3333333 3.6666667
2 0.8279142 3.3333333
3 -0.3280125 0.8279142
4 -0.3280125 -0.3280125
5 -4.3248593 -0.3280125
6 -0.6926295 -4.3248593
7 -0.1442503 -0.6926295
8 0.6826292 -0.1442503
9 1.6213424 0.6826292
10 1.7837229 1.6213424
11 4.2417975 1.7837229
12 1.5374917 4.2417975
13 2.2041584 1.5374917
14 0.6665192 2.2041584
15 -1.7474625 0.6665192
16 -1.7474625 -1.7474625
17 -0.8087493 -1.7474625
18 1.3073705 -0.8087493
19 2.8557497 1.3073705
20 3.6826292 2.8557497
21 2.6213424 3.6826292
22 2.7837229 2.6213424
23 2.2417975 2.7837229
24 -4.4625083 2.2417975
25 -3.7958416 -4.4625083
26 -3.3334808 -3.7958416
27 -4.7474625 -3.3334808
28 -3.7474625 -4.7474625
29 -3.8087493 -3.7474625
30 -2.6926295 -3.8087493
31 -2.1442503 -2.6926295
32 -0.3173708 -2.1442503
33 0.6213424 -0.3173708
34 0.7837229 0.6213424
35 2.3579665 0.7837229
36 -5.7527833 2.3579665
37 -4.0861166 -5.7527833
38 -1.6559758 -4.0861166
39 -0.3280125 -1.6559758
40 2.6719875 -0.3280125
41 4.6881467 2.6719875
42 3.4365455 4.6881467
43 4.6559267 3.4365455
44 6.5215292 4.6559267
45 4.9828553 6.5215292
46 5.3323478 4.9828553
47 3.8678685 5.3323478
48 4.7958416 3.8678685
49 2.4625083 4.7958416
50 4.1507041 2.4625083
51 6.2206124 4.1507041
52 5.2206124 6.2206124
53 4.3851606 5.2206124
54 1.9851704 4.3851606
55 0.5335496 1.9851704
56 -3.6395709 0.5335496
57 -1.7008577 -3.6395709
58 -0.5384772 -1.7008577
59 1.9195974 -0.5384772
60 0.2152916 1.9195974
61 -0.1180417 0.2152916
62 -0.6556809 -0.1180417
63 0.9303374 -0.6556809
64 -2.0696626 0.9303374
65 -0.1309494 -2.0696626
66 -3.3438276 -0.1309494
67 -5.7567254 -3.3438276
68 -6.9298459 -5.7567254
69 -8.1460247 -6.9298459
70 -10.1450392 -8.1460247
71 -14.6290275 -10.1450392
> 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/70jdm1258711054.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/8z5h11258711054.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/9xgs41258711054.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/106imq1258711054.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/111alv1258711054.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/12eeaq1258711054.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/136alq1258711054.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/14c7i61258711054.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/15hojo1258711054.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/16p0rk1258711054.tab")
+ }
>
> system("convert tmp/1bysl1258711054.ps tmp/1bysl1258711054.png")
> system("convert tmp/2xj5c1258711054.ps tmp/2xj5c1258711054.png")
> system("convert tmp/3900q1258711054.ps tmp/3900q1258711054.png")
> system("convert tmp/4a9jj1258711054.ps tmp/4a9jj1258711054.png")
> system("convert tmp/5h2s81258711054.ps tmp/5h2s81258711054.png")
> system("convert tmp/627bg1258711054.ps tmp/627bg1258711054.png")
> system("convert tmp/70jdm1258711054.ps tmp/70jdm1258711054.png")
> system("convert tmp/8z5h11258711054.ps tmp/8z5h11258711054.png")
> system("convert tmp/9xgs41258711054.ps tmp/9xgs41258711054.png")
> system("convert tmp/106imq1258711054.ps tmp/106imq1258711054.png")
>
>
> proc.time()
user system elapsed
2.523 1.558 3.381