R version 2.8.0 (2008-10-20)
Copyright (C) 2008 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(95.2,0,95.00,0,94.00,0,92.2,0,91.00,0,91.2,0,103.4,1,105.00,1,104.6,1,103.8,0,101.8,0,102.4,0,103.8,0,103.4,0,102.00,0,101.8,0,100.2,0,101.4,0,113.8,1,116.00,1,115.6,1,113.00,0,109.4,0,111.00,0,112.4,0,112.2,0,111.00,0,108.8,0,107.4,0,108.6,0,118.8,1,122.2,1,122.6,1,122.2,0,118.8,0,119.00,0,118.2,0,117.8,0,116.8,0,114.6,0,113.4,0,113.8,0,124.2,1,125.8,1,125.6,1,122.4,0,119.00,0,119.4,0,118.6,0,118.00,0,116.00,0,114.8,0,114.6,0,114.6,0,124.00,1,125.2,1,124.00,1,117.6,0,113.2,0,111.4,0,112.2,0,109.8,0,106.4,0,105.2,0,102.2,0,99.8,0,111.00,1,113.00,1,108.4,1,105.4,0,102.00,0,102.8,0),dim=c(2,72),dimnames=list(c('Werkloosheid','dummy'),1:72))
> y <- array(NA,dim=c(2,72),dimnames=list(c('Werkloosheid','dummy'),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 = 'Do not include Seasonal 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 dummy
1 95.2 0
2 95.0 0
3 94.0 0
4 92.2 0
5 91.0 0
6 91.2 0
7 103.4 1
8 105.0 1
9 104.6 1
10 103.8 0
11 101.8 0
12 102.4 0
13 103.8 0
14 103.4 0
15 102.0 0
16 101.8 0
17 100.2 0
18 101.4 0
19 113.8 1
20 116.0 1
21 115.6 1
22 113.0 0
23 109.4 0
24 111.0 0
25 112.4 0
26 112.2 0
27 111.0 0
28 108.8 0
29 107.4 0
30 108.6 0
31 118.8 1
32 122.2 1
33 122.6 1
34 122.2 0
35 118.8 0
36 119.0 0
37 118.2 0
38 117.8 0
39 116.8 0
40 114.6 0
41 113.4 0
42 113.8 0
43 124.2 1
44 125.8 1
45 125.6 1
46 122.4 0
47 119.0 0
48 119.4 0
49 118.6 0
50 118.0 0
51 116.0 0
52 114.8 0
53 114.6 0
54 114.6 0
55 124.0 1
56 125.2 1
57 124.0 1
58 117.6 0
59 113.2 0
60 111.4 0
61 112.2 0
62 109.8 0
63 106.4 0
64 105.2 0
65 102.2 0
66 99.8 0
67 111.0 1
68 113.0 1
69 108.4 1
70 105.4 0
71 102.0 0
72 102.8 0
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) dummy
108.759 8.085
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-17.759 -6.409 1.498 7.156 13.641
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 108.759 1.125 96.691 < 2e-16 ***
dummy 8.085 2.250 3.594 0.000602 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 8.266 on 70 degrees of freedom
Multiple R-squared: 0.1558, Adjusted R-squared: 0.1437
F-statistic: 12.92 on 1 and 70 DF, p-value: 0.0006016
> 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.0308891149 0.0617782298 0.9691108851
[2,] 0.0149920629 0.0299841257 0.9850079371
[3,] 0.0043745436 0.0087490872 0.9956254564
[4,] 0.0014997484 0.0029994967 0.9985002516
[5,] 0.0004614344 0.0009228687 0.9995385656
[6,] 0.0588203660 0.1176407321 0.9411796340
[7,] 0.0927639161 0.1855278322 0.9072360839
[8,] 0.1218325916 0.2436651831 0.8781674084
[9,] 0.1628070220 0.3256140439 0.8371929780
[10,] 0.1798220399 0.3596440798 0.8201779601
[11,] 0.1712537935 0.3425075869 0.8287462065
[12,] 0.1616521114 0.3233042227 0.8383478886
[13,] 0.1516783019 0.3033566038 0.8483216981
[14,] 0.1490118413 0.2980236827 0.8509881587
[15,] 0.2153417577 0.4306835154 0.7846582423
[16,] 0.2773512273 0.5547024546 0.7226487727
[17,] 0.2905870131 0.5811740262 0.7094129869
[18,] 0.5319591068 0.9360817864 0.4680408932
[19,] 0.5929029805 0.8141940390 0.4070970195
[20,] 0.6599939981 0.6800120038 0.3400060019
[21,] 0.7251707441 0.5496585118 0.2748292559
[22,] 0.7599941753 0.4800116495 0.2400058247
[23,] 0.7648006839 0.4703986321 0.2351993161
[24,] 0.7470457992 0.5059084017 0.2529542008
[25,] 0.7239083965 0.5521832070 0.2760916035
[26,] 0.7012074488 0.5975851024 0.2987925512
[27,] 0.6937396636 0.6125206729 0.3062603364
[28,] 0.7103450498 0.5793099004 0.2896549502
[29,] 0.7128211298 0.5743577403 0.2871788702
[30,] 0.8679237630 0.2641524739 0.1320762370
[31,] 0.9052699496 0.1894601009 0.0947300504
[32,] 0.9299875479 0.1400249042 0.0700124521
[33,] 0.9409735898 0.1180528203 0.0590264102
[34,] 0.9462807953 0.1074384095 0.0537192047
[35,] 0.9449917860 0.1100164279 0.0550082140
[36,] 0.9335408130 0.1329183741 0.0664591870
[37,] 0.9153035634 0.1693928733 0.0846964366
[38,] 0.8943131680 0.2113736640 0.1056868320
[39,] 0.8850087256 0.2299825488 0.1149912744
[40,] 0.8866077020 0.2267845960 0.1133922980
[41,] 0.8885860617 0.2228278767 0.1114139383
[42,] 0.9316118708 0.1367762584 0.0683881292
[43,] 0.9393781454 0.1212437092 0.0606218546
[44,] 0.9503374047 0.0993251905 0.0496625953
[45,] 0.9572353393 0.0855293214 0.0427646607
[46,] 0.9624536249 0.0750927503 0.0375463751
[47,] 0.9601534089 0.0796931821 0.0398465911
[48,] 0.9539147632 0.0921704736 0.0460852368
[49,] 0.9477265497 0.1045469005 0.0522734503
[50,] 0.9435651761 0.1128696479 0.0564348239
[51,] 0.9408945435 0.1182109131 0.0591054565
[52,] 0.9593249214 0.0813501573 0.0406750786
[53,] 0.9832140647 0.0335718706 0.0167859353
[54,] 0.9952281792 0.0095436416 0.0047718208
[55,] 0.9966820736 0.0066358528 0.0033179264
[56,] 0.9970328989 0.0059342023 0.0029671011
[57,] 0.9990967021 0.0018065958 0.0009032979
[58,] 0.9996917526 0.0006164947 0.0003082474
[59,] 0.9995258479 0.0009483042 0.0004741521
[60,] 0.9989514315 0.0020971370 0.0010485685
[61,] 0.9957317235 0.0085365530 0.0042682765
[62,] 0.9932338551 0.0135322898 0.0067661449
[63,] 0.9707586383 0.0584827235 0.0292413617
> postscript(file="/var/www/html/rcomp/tmp/1r65b1228656114.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/2u2ka1228656114.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/39q2s1228656114.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/4fnvs1228656114.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/5zhid1228656114.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
-13.55925926 -13.75925926 -14.75925926 -16.55925926 -17.75925926 -17.55925926
7 8 9 10 11 12
-13.44444444 -11.84444444 -12.24444444 -4.95925926 -6.95925926 -6.35925926
13 14 15 16 17 18
-4.95925926 -5.35925926 -6.75925926 -6.95925926 -8.55925926 -7.35925926
19 20 21 22 23 24
-3.04444444 -0.84444444 -1.24444444 4.24074074 0.64074074 2.24074074
25 26 27 28 29 30
3.64074074 3.44074074 2.24074074 0.04074074 -1.35925926 -0.15925926
31 32 33 34 35 36
1.95555556 5.35555556 5.75555556 13.44074074 10.04074074 10.24074074
37 38 39 40 41 42
9.44074074 9.04074074 8.04074074 5.84074074 4.64074074 5.04074074
43 44 45 46 47 48
7.35555556 8.95555556 8.75555556 13.64074074 10.24074074 10.64074074
49 50 51 52 53 54
9.84074074 9.24074074 7.24074074 6.04074074 5.84074074 5.84074074
55 56 57 58 59 60
7.15555556 8.35555556 7.15555556 8.84074074 4.44074074 2.64074074
61 62 63 64 65 66
3.44074074 1.04074074 -2.35925926 -3.55925926 -6.55925926 -8.95925926
67 68 69 70 71 72
-5.84444444 -3.84444444 -8.44444444 -3.35925926 -6.75925926 -5.95925926
> postscript(file="/var/www/html/rcomp/tmp/6mtt11228656114.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 -13.55925926 NA
1 -13.75925926 -13.55925926
2 -14.75925926 -13.75925926
3 -16.55925926 -14.75925926
4 -17.75925926 -16.55925926
5 -17.55925926 -17.75925926
6 -13.44444444 -17.55925926
7 -11.84444444 -13.44444444
8 -12.24444444 -11.84444444
9 -4.95925926 -12.24444444
10 -6.95925926 -4.95925926
11 -6.35925926 -6.95925926
12 -4.95925926 -6.35925926
13 -5.35925926 -4.95925926
14 -6.75925926 -5.35925926
15 -6.95925926 -6.75925926
16 -8.55925926 -6.95925926
17 -7.35925926 -8.55925926
18 -3.04444444 -7.35925926
19 -0.84444444 -3.04444444
20 -1.24444444 -0.84444444
21 4.24074074 -1.24444444
22 0.64074074 4.24074074
23 2.24074074 0.64074074
24 3.64074074 2.24074074
25 3.44074074 3.64074074
26 2.24074074 3.44074074
27 0.04074074 2.24074074
28 -1.35925926 0.04074074
29 -0.15925926 -1.35925926
30 1.95555556 -0.15925926
31 5.35555556 1.95555556
32 5.75555556 5.35555556
33 13.44074074 5.75555556
34 10.04074074 13.44074074
35 10.24074074 10.04074074
36 9.44074074 10.24074074
37 9.04074074 9.44074074
38 8.04074074 9.04074074
39 5.84074074 8.04074074
40 4.64074074 5.84074074
41 5.04074074 4.64074074
42 7.35555556 5.04074074
43 8.95555556 7.35555556
44 8.75555556 8.95555556
45 13.64074074 8.75555556
46 10.24074074 13.64074074
47 10.64074074 10.24074074
48 9.84074074 10.64074074
49 9.24074074 9.84074074
50 7.24074074 9.24074074
51 6.04074074 7.24074074
52 5.84074074 6.04074074
53 5.84074074 5.84074074
54 7.15555556 5.84074074
55 8.35555556 7.15555556
56 7.15555556 8.35555556
57 8.84074074 7.15555556
58 4.44074074 8.84074074
59 2.64074074 4.44074074
60 3.44074074 2.64074074
61 1.04074074 3.44074074
62 -2.35925926 1.04074074
63 -3.55925926 -2.35925926
64 -6.55925926 -3.55925926
65 -8.95925926 -6.55925926
66 -5.84444444 -8.95925926
67 -3.84444444 -5.84444444
68 -8.44444444 -3.84444444
69 -3.35925926 -8.44444444
70 -6.75925926 -3.35925926
71 -5.95925926 -6.75925926
72 NA -5.95925926
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -13.75925926 -13.55925926
[2,] -14.75925926 -13.75925926
[3,] -16.55925926 -14.75925926
[4,] -17.75925926 -16.55925926
[5,] -17.55925926 -17.75925926
[6,] -13.44444444 -17.55925926
[7,] -11.84444444 -13.44444444
[8,] -12.24444444 -11.84444444
[9,] -4.95925926 -12.24444444
[10,] -6.95925926 -4.95925926
[11,] -6.35925926 -6.95925926
[12,] -4.95925926 -6.35925926
[13,] -5.35925926 -4.95925926
[14,] -6.75925926 -5.35925926
[15,] -6.95925926 -6.75925926
[16,] -8.55925926 -6.95925926
[17,] -7.35925926 -8.55925926
[18,] -3.04444444 -7.35925926
[19,] -0.84444444 -3.04444444
[20,] -1.24444444 -0.84444444
[21,] 4.24074074 -1.24444444
[22,] 0.64074074 4.24074074
[23,] 2.24074074 0.64074074
[24,] 3.64074074 2.24074074
[25,] 3.44074074 3.64074074
[26,] 2.24074074 3.44074074
[27,] 0.04074074 2.24074074
[28,] -1.35925926 0.04074074
[29,] -0.15925926 -1.35925926
[30,] 1.95555556 -0.15925926
[31,] 5.35555556 1.95555556
[32,] 5.75555556 5.35555556
[33,] 13.44074074 5.75555556
[34,] 10.04074074 13.44074074
[35,] 10.24074074 10.04074074
[36,] 9.44074074 10.24074074
[37,] 9.04074074 9.44074074
[38,] 8.04074074 9.04074074
[39,] 5.84074074 8.04074074
[40,] 4.64074074 5.84074074
[41,] 5.04074074 4.64074074
[42,] 7.35555556 5.04074074
[43,] 8.95555556 7.35555556
[44,] 8.75555556 8.95555556
[45,] 13.64074074 8.75555556
[46,] 10.24074074 13.64074074
[47,] 10.64074074 10.24074074
[48,] 9.84074074 10.64074074
[49,] 9.24074074 9.84074074
[50,] 7.24074074 9.24074074
[51,] 6.04074074 7.24074074
[52,] 5.84074074 6.04074074
[53,] 5.84074074 5.84074074
[54,] 7.15555556 5.84074074
[55,] 8.35555556 7.15555556
[56,] 7.15555556 8.35555556
[57,] 8.84074074 7.15555556
[58,] 4.44074074 8.84074074
[59,] 2.64074074 4.44074074
[60,] 3.44074074 2.64074074
[61,] 1.04074074 3.44074074
[62,] -2.35925926 1.04074074
[63,] -3.55925926 -2.35925926
[64,] -6.55925926 -3.55925926
[65,] -8.95925926 -6.55925926
[66,] -5.84444444 -8.95925926
[67,] -3.84444444 -5.84444444
[68,] -8.44444444 -3.84444444
[69,] -3.35925926 -8.44444444
[70,] -6.75925926 -3.35925926
[71,] -5.95925926 -6.75925926
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -13.75925926 -13.55925926
2 -14.75925926 -13.75925926
3 -16.55925926 -14.75925926
4 -17.75925926 -16.55925926
5 -17.55925926 -17.75925926
6 -13.44444444 -17.55925926
7 -11.84444444 -13.44444444
8 -12.24444444 -11.84444444
9 -4.95925926 -12.24444444
10 -6.95925926 -4.95925926
11 -6.35925926 -6.95925926
12 -4.95925926 -6.35925926
13 -5.35925926 -4.95925926
14 -6.75925926 -5.35925926
15 -6.95925926 -6.75925926
16 -8.55925926 -6.95925926
17 -7.35925926 -8.55925926
18 -3.04444444 -7.35925926
19 -0.84444444 -3.04444444
20 -1.24444444 -0.84444444
21 4.24074074 -1.24444444
22 0.64074074 4.24074074
23 2.24074074 0.64074074
24 3.64074074 2.24074074
25 3.44074074 3.64074074
26 2.24074074 3.44074074
27 0.04074074 2.24074074
28 -1.35925926 0.04074074
29 -0.15925926 -1.35925926
30 1.95555556 -0.15925926
31 5.35555556 1.95555556
32 5.75555556 5.35555556
33 13.44074074 5.75555556
34 10.04074074 13.44074074
35 10.24074074 10.04074074
36 9.44074074 10.24074074
37 9.04074074 9.44074074
38 8.04074074 9.04074074
39 5.84074074 8.04074074
40 4.64074074 5.84074074
41 5.04074074 4.64074074
42 7.35555556 5.04074074
43 8.95555556 7.35555556
44 8.75555556 8.95555556
45 13.64074074 8.75555556
46 10.24074074 13.64074074
47 10.64074074 10.24074074
48 9.84074074 10.64074074
49 9.24074074 9.84074074
50 7.24074074 9.24074074
51 6.04074074 7.24074074
52 5.84074074 6.04074074
53 5.84074074 5.84074074
54 7.15555556 5.84074074
55 8.35555556 7.15555556
56 7.15555556 8.35555556
57 8.84074074 7.15555556
58 4.44074074 8.84074074
59 2.64074074 4.44074074
60 3.44074074 2.64074074
61 1.04074074 3.44074074
62 -2.35925926 1.04074074
63 -3.55925926 -2.35925926
64 -6.55925926 -3.55925926
65 -8.95925926 -6.55925926
66 -5.84444444 -8.95925926
67 -3.84444444 -5.84444444
68 -8.44444444 -3.84444444
69 -3.35925926 -8.44444444
70 -6.75925926 -3.35925926
71 -5.95925926 -6.75925926
> 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/7vdlm1228656114.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/8g3rw1228656114.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/97ogq1228656114.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/10lh111228656114.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/11c43h1228656114.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/12ehqq1228656114.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/13qd5g1228656114.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/14tu8k1228656114.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/1512al1228656114.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/16e90t1228656114.tab")
+ }
>
> system("convert tmp/1r65b1228656114.ps tmp/1r65b1228656114.png")
> system("convert tmp/2u2ka1228656114.ps tmp/2u2ka1228656114.png")
> system("convert tmp/39q2s1228656114.ps tmp/39q2s1228656114.png")
> system("convert tmp/4fnvs1228656114.ps tmp/4fnvs1228656114.png")
> system("convert tmp/5zhid1228656114.ps tmp/5zhid1228656114.png")
> system("convert tmp/6mtt11228656114.ps tmp/6mtt11228656114.png")
> system("convert tmp/7vdlm1228656114.ps tmp/7vdlm1228656114.png")
> system("convert tmp/8g3rw1228656114.ps tmp/8g3rw1228656114.png")
> system("convert tmp/97ogq1228656114.ps tmp/97ogq1228656114.png")
> system("convert tmp/10lh111228656114.ps tmp/10lh111228656114.png")
>
>
> proc.time()
user system elapsed
2.614 1.602 3.263