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(269645,0,267037,0,258113,0,262813,0,267413,0,267366,0,264777,0,258863,0,254844,0,254868,0,277267,0,285351,0,286602,0,283042,0,276687,0,277915,0,277128,0,277103,0,275037,0,270150,0,267140,0,264993,0,287259,0,291186,0,292300,0,288186,0,281477,0,282656,1,280190,1,280408,1,276836,1,275216,1,274352,1,271311,1,289802,1,290726,1,292300,1,278506,1,269826,1,265861,1,269034,1,264176,1,255198,1,253353,1,246057,1,235372,1,258556,1,260993,1,254663,1,250643,1,243422,1,247105,1,248541,1,245039,1,237080,1,237085,1,225554,1,226839,1,247934,1,248333,1,246969,1,245098,1,246263,1),dim=c(2,63),dimnames=list(c('Mannen','Dummy'),1:63))
> y <- array(NA,dim=c(2,63),dimnames=list(c('Mannen','Dummy'),1:63))
> 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
Mannen Dummy M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 269645 0 1 0 0 0 0 0 0 0 0 0 0 1
2 267037 0 0 1 0 0 0 0 0 0 0 0 0 2
3 258113 0 0 0 1 0 0 0 0 0 0 0 0 3
4 262813 0 0 0 0 1 0 0 0 0 0 0 0 4
5 267413 0 0 0 0 0 1 0 0 0 0 0 0 5
6 267366 0 0 0 0 0 0 1 0 0 0 0 0 6
7 264777 0 0 0 0 0 0 0 1 0 0 0 0 7
8 258863 0 0 0 0 0 0 0 0 1 0 0 0 8
9 254844 0 0 0 0 0 0 0 0 0 1 0 0 9
10 254868 0 0 0 0 0 0 0 0 0 0 1 0 10
11 277267 0 0 0 0 0 0 0 0 0 0 0 1 11
12 285351 0 0 0 0 0 0 0 0 0 0 0 0 12
13 286602 0 1 0 0 0 0 0 0 0 0 0 0 13
14 283042 0 0 1 0 0 0 0 0 0 0 0 0 14
15 276687 0 0 0 1 0 0 0 0 0 0 0 0 15
16 277915 0 0 0 0 1 0 0 0 0 0 0 0 16
17 277128 0 0 0 0 0 1 0 0 0 0 0 0 17
18 277103 0 0 0 0 0 0 1 0 0 0 0 0 18
19 275037 0 0 0 0 0 0 0 1 0 0 0 0 19
20 270150 0 0 0 0 0 0 0 0 1 0 0 0 20
21 267140 0 0 0 0 0 0 0 0 0 1 0 0 21
22 264993 0 0 0 0 0 0 0 0 0 0 1 0 22
23 287259 0 0 0 0 0 0 0 0 0 0 0 1 23
24 291186 0 0 0 0 0 0 0 0 0 0 0 0 24
25 292300 0 1 0 0 0 0 0 0 0 0 0 0 25
26 288186 0 0 1 0 0 0 0 0 0 0 0 0 26
27 281477 0 0 0 1 0 0 0 0 0 0 0 0 27
28 282656 1 0 0 0 1 0 0 0 0 0 0 0 28
29 280190 1 0 0 0 0 1 0 0 0 0 0 0 29
30 280408 1 0 0 0 0 0 1 0 0 0 0 0 30
31 276836 1 0 0 0 0 0 0 1 0 0 0 0 31
32 275216 1 0 0 0 0 0 0 0 1 0 0 0 32
33 274352 1 0 0 0 0 0 0 0 0 1 0 0 33
34 271311 1 0 0 0 0 0 0 0 0 0 1 0 34
35 289802 1 0 0 0 0 0 0 0 0 0 0 1 35
36 290726 1 0 0 0 0 0 0 0 0 0 0 0 36
37 292300 1 1 0 0 0 0 0 0 0 0 0 0 37
38 278506 1 0 1 0 0 0 0 0 0 0 0 0 38
39 269826 1 0 0 1 0 0 0 0 0 0 0 0 39
40 265861 1 0 0 0 1 0 0 0 0 0 0 0 40
41 269034 1 0 0 0 0 1 0 0 0 0 0 0 41
42 264176 1 0 0 0 0 0 1 0 0 0 0 0 42
43 255198 1 0 0 0 0 0 0 1 0 0 0 0 43
44 253353 1 0 0 0 0 0 0 0 1 0 0 0 44
45 246057 1 0 0 0 0 0 0 0 0 1 0 0 45
46 235372 1 0 0 0 0 0 0 0 0 0 1 0 46
47 258556 1 0 0 0 0 0 0 0 0 0 0 1 47
48 260993 1 0 0 0 0 0 0 0 0 0 0 0 48
49 254663 1 1 0 0 0 0 0 0 0 0 0 0 49
50 250643 1 0 1 0 0 0 0 0 0 0 0 0 50
51 243422 1 0 0 1 0 0 0 0 0 0 0 0 51
52 247105 1 0 0 0 1 0 0 0 0 0 0 0 52
53 248541 1 0 0 0 0 1 0 0 0 0 0 0 53
54 245039 1 0 0 0 0 0 1 0 0 0 0 0 54
55 237080 1 0 0 0 0 0 0 1 0 0 0 0 55
56 237085 1 0 0 0 0 0 0 0 1 0 0 0 56
57 225554 1 0 0 0 0 0 0 0 0 1 0 0 57
58 226839 1 0 0 0 0 0 0 0 0 0 1 0 58
59 247934 1 0 0 0 0 0 0 0 0 0 0 1 59
60 248333 1 0 0 0 0 0 0 0 0 0 0 0 60
61 246969 1 1 0 0 0 0 0 0 0 0 0 0 61
62 245098 1 0 1 0 0 0 0 0 0 0 0 0 62
63 246263 1 0 0 1 0 0 0 0 0 0 0 0 63
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Dummy M1 M2 M3 M4
298214.3 12494.8 -4543.1 -8693.4 -13969.8 -14801.9
M5 M6 M7 M8 M9 M10
-12766.4 -13565.0 -17753.5 -19761.4 -24261.2 -26329.7
M11 t
-3998.5 -844.3
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-23598.8 -9175.6 -447.2 9322.1 20615.8
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 298214.3 6725.5 44.341 < 2e-16 ***
Dummy 12494.8 6529.2 1.914 0.06151 .
M1 -4543.1 7710.0 -0.589 0.55840
M2 -8693.4 7704.6 -1.128 0.26467
M3 -13969.8 7703.3 -1.813 0.07588 .
M4 -14801.9 8164.1 -1.813 0.07595 .
M5 -12766.4 8134.9 -1.569 0.12300
M6 -13565.0 8109.4 -1.673 0.10075
M7 -17753.5 8087.8 -2.195 0.03292 *
M8 -19761.4 8070.1 -2.449 0.01796 *
M9 -24261.2 8056.3 -3.011 0.00410 **
M10 -26329.7 8046.5 -3.272 0.00196 **
M11 -3998.5 8040.5 -0.497 0.62121
t -844.3 178.3 -4.735 1.90e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 12710 on 49 degrees of freedom
Multiple R-squared: 0.5616, Adjusted R-squared: 0.4453
F-statistic: 4.828 on 13 and 49 DF, p-value: 2.508e-05
> 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.131814e-01 2.263627e-01 0.88681863
[2,] 9.486807e-02 1.897361e-01 0.90513193
[3,] 5.817207e-02 1.163441e-01 0.94182793
[4,] 3.180961e-02 6.361921e-02 0.96819039
[5,] 1.437071e-02 2.874143e-02 0.98562929
[6,] 8.039208e-03 1.607842e-02 0.99196079
[7,] 4.338308e-03 8.676616e-03 0.99566169
[8,] 6.462000e-03 1.292400e-02 0.99353800
[9,] 4.200207e-03 8.400415e-03 0.99579979
[10,] 2.869476e-03 5.738953e-03 0.99713052
[11,] 1.407018e-03 2.814036e-03 0.99859298
[12,] 5.562654e-04 1.112531e-03 0.99944373
[13,] 3.835484e-04 7.670967e-04 0.99961645
[14,] 1.647873e-04 3.295746e-04 0.99983521
[15,] 7.026652e-05 1.405330e-04 0.99992973
[16,] 2.650884e-05 5.301768e-05 0.99997349
[17,] 2.650832e-05 5.301664e-05 0.99997349
[18,] 2.523380e-05 5.046760e-05 0.99997477
[19,] 2.187737e-05 4.375473e-05 0.99997812
[20,] 1.087604e-04 2.175209e-04 0.99989124
[21,] 2.195802e-03 4.391603e-03 0.99780420
[22,] 6.302398e-02 1.260480e-01 0.93697602
[23,] 1.712064e-01 3.424127e-01 0.82879364
[24,] 4.798351e-01 9.596701e-01 0.52016493
[25,] 6.215068e-01 7.569864e-01 0.37849320
[26,] 7.374924e-01 5.250151e-01 0.26250757
[27,] 8.242673e-01 3.514654e-01 0.17573269
[28,] 8.258049e-01 3.483903e-01 0.17419514
[29,] 9.238493e-01 1.523014e-01 0.07615072
[30,] 8.726327e-01 2.547346e-01 0.12736730
> postscript(file="/var/www/html/rcomp/tmp/1xqr91229461290.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/26w7e1229461290.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/3bbdx1229461290.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/4bo3o1229461290.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/5p90s1229461290.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 = 63
Frequency = 1
1 2 3 4 5 6
-23181.9444 -20795.4444 -23598.7778 -17222.3917 -13813.5917 -12217.7917
7 8 9 10 11 12
-9773.9917 -12835.7917 -11510.7917 -8573.9917 -7661.9917 -2732.1917
13 14 15 16 17 18
3906.1972 5340.6972 5106.3639 8010.7500 6032.5500 7650.3500
19 20 21 22 23 24
10617.1500 8582.3500 10916.3500 11682.1500 12461.1500 13233.9500
25 26 27 28 29 30
19735.3389 20615.8389 20027.5056 10388.0722 6730.8722 8591.6722
31 32 33 34 35 36
10052.4722 11284.6722 15764.6722 15636.4722 12640.4722 10410.2722
37 38 39 40 41 42
17371.6611 8572.1611 6012.8278 3724.2139 5706.0139 2490.8139
43 44 45 46 47 48
-1454.3861 -447.1861 -2399.1861 -10171.3861 -8474.3861 -9191.5861
49 50 51 52 53 54
-10134.1972 -9159.6972 -10260.0306 -4900.6444 -4655.8444 -6515.0444
55 56 57 58 59 60
-9441.2444 -6584.0444 -12771.0444 -8573.2444 -8965.2444 -11720.4444
61 62 63
-7697.0556 -4573.5556 2712.1111
> postscript(file="/var/www/html/rcomp/tmp/6ro2f1229461290.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 = 63
Frequency = 1
lag(myerror, k = 1) myerror
0 -23181.9444 NA
1 -20795.4444 -23181.9444
2 -23598.7778 -20795.4444
3 -17222.3917 -23598.7778
4 -13813.5917 -17222.3917
5 -12217.7917 -13813.5917
6 -9773.9917 -12217.7917
7 -12835.7917 -9773.9917
8 -11510.7917 -12835.7917
9 -8573.9917 -11510.7917
10 -7661.9917 -8573.9917
11 -2732.1917 -7661.9917
12 3906.1972 -2732.1917
13 5340.6972 3906.1972
14 5106.3639 5340.6972
15 8010.7500 5106.3639
16 6032.5500 8010.7500
17 7650.3500 6032.5500
18 10617.1500 7650.3500
19 8582.3500 10617.1500
20 10916.3500 8582.3500
21 11682.1500 10916.3500
22 12461.1500 11682.1500
23 13233.9500 12461.1500
24 19735.3389 13233.9500
25 20615.8389 19735.3389
26 20027.5056 20615.8389
27 10388.0722 20027.5056
28 6730.8722 10388.0722
29 8591.6722 6730.8722
30 10052.4722 8591.6722
31 11284.6722 10052.4722
32 15764.6722 11284.6722
33 15636.4722 15764.6722
34 12640.4722 15636.4722
35 10410.2722 12640.4722
36 17371.6611 10410.2722
37 8572.1611 17371.6611
38 6012.8278 8572.1611
39 3724.2139 6012.8278
40 5706.0139 3724.2139
41 2490.8139 5706.0139
42 -1454.3861 2490.8139
43 -447.1861 -1454.3861
44 -2399.1861 -447.1861
45 -10171.3861 -2399.1861
46 -8474.3861 -10171.3861
47 -9191.5861 -8474.3861
48 -10134.1972 -9191.5861
49 -9159.6972 -10134.1972
50 -10260.0306 -9159.6972
51 -4900.6444 -10260.0306
52 -4655.8444 -4900.6444
53 -6515.0444 -4655.8444
54 -9441.2444 -6515.0444
55 -6584.0444 -9441.2444
56 -12771.0444 -6584.0444
57 -8573.2444 -12771.0444
58 -8965.2444 -8573.2444
59 -11720.4444 -8965.2444
60 -7697.0556 -11720.4444
61 -4573.5556 -7697.0556
62 2712.1111 -4573.5556
63 NA 2712.1111
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -20795.4444 -23181.9444
[2,] -23598.7778 -20795.4444
[3,] -17222.3917 -23598.7778
[4,] -13813.5917 -17222.3917
[5,] -12217.7917 -13813.5917
[6,] -9773.9917 -12217.7917
[7,] -12835.7917 -9773.9917
[8,] -11510.7917 -12835.7917
[9,] -8573.9917 -11510.7917
[10,] -7661.9917 -8573.9917
[11,] -2732.1917 -7661.9917
[12,] 3906.1972 -2732.1917
[13,] 5340.6972 3906.1972
[14,] 5106.3639 5340.6972
[15,] 8010.7500 5106.3639
[16,] 6032.5500 8010.7500
[17,] 7650.3500 6032.5500
[18,] 10617.1500 7650.3500
[19,] 8582.3500 10617.1500
[20,] 10916.3500 8582.3500
[21,] 11682.1500 10916.3500
[22,] 12461.1500 11682.1500
[23,] 13233.9500 12461.1500
[24,] 19735.3389 13233.9500
[25,] 20615.8389 19735.3389
[26,] 20027.5056 20615.8389
[27,] 10388.0722 20027.5056
[28,] 6730.8722 10388.0722
[29,] 8591.6722 6730.8722
[30,] 10052.4722 8591.6722
[31,] 11284.6722 10052.4722
[32,] 15764.6722 11284.6722
[33,] 15636.4722 15764.6722
[34,] 12640.4722 15636.4722
[35,] 10410.2722 12640.4722
[36,] 17371.6611 10410.2722
[37,] 8572.1611 17371.6611
[38,] 6012.8278 8572.1611
[39,] 3724.2139 6012.8278
[40,] 5706.0139 3724.2139
[41,] 2490.8139 5706.0139
[42,] -1454.3861 2490.8139
[43,] -447.1861 -1454.3861
[44,] -2399.1861 -447.1861
[45,] -10171.3861 -2399.1861
[46,] -8474.3861 -10171.3861
[47,] -9191.5861 -8474.3861
[48,] -10134.1972 -9191.5861
[49,] -9159.6972 -10134.1972
[50,] -10260.0306 -9159.6972
[51,] -4900.6444 -10260.0306
[52,] -4655.8444 -4900.6444
[53,] -6515.0444 -4655.8444
[54,] -9441.2444 -6515.0444
[55,] -6584.0444 -9441.2444
[56,] -12771.0444 -6584.0444
[57,] -8573.2444 -12771.0444
[58,] -8965.2444 -8573.2444
[59,] -11720.4444 -8965.2444
[60,] -7697.0556 -11720.4444
[61,] -4573.5556 -7697.0556
[62,] 2712.1111 -4573.5556
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -20795.4444 -23181.9444
2 -23598.7778 -20795.4444
3 -17222.3917 -23598.7778
4 -13813.5917 -17222.3917
5 -12217.7917 -13813.5917
6 -9773.9917 -12217.7917
7 -12835.7917 -9773.9917
8 -11510.7917 -12835.7917
9 -8573.9917 -11510.7917
10 -7661.9917 -8573.9917
11 -2732.1917 -7661.9917
12 3906.1972 -2732.1917
13 5340.6972 3906.1972
14 5106.3639 5340.6972
15 8010.7500 5106.3639
16 6032.5500 8010.7500
17 7650.3500 6032.5500
18 10617.1500 7650.3500
19 8582.3500 10617.1500
20 10916.3500 8582.3500
21 11682.1500 10916.3500
22 12461.1500 11682.1500
23 13233.9500 12461.1500
24 19735.3389 13233.9500
25 20615.8389 19735.3389
26 20027.5056 20615.8389
27 10388.0722 20027.5056
28 6730.8722 10388.0722
29 8591.6722 6730.8722
30 10052.4722 8591.6722
31 11284.6722 10052.4722
32 15764.6722 11284.6722
33 15636.4722 15764.6722
34 12640.4722 15636.4722
35 10410.2722 12640.4722
36 17371.6611 10410.2722
37 8572.1611 17371.6611
38 6012.8278 8572.1611
39 3724.2139 6012.8278
40 5706.0139 3724.2139
41 2490.8139 5706.0139
42 -1454.3861 2490.8139
43 -447.1861 -1454.3861
44 -2399.1861 -447.1861
45 -10171.3861 -2399.1861
46 -8474.3861 -10171.3861
47 -9191.5861 -8474.3861
48 -10134.1972 -9191.5861
49 -9159.6972 -10134.1972
50 -10260.0306 -9159.6972
51 -4900.6444 -10260.0306
52 -4655.8444 -4900.6444
53 -6515.0444 -4655.8444
54 -9441.2444 -6515.0444
55 -6584.0444 -9441.2444
56 -12771.0444 -6584.0444
57 -8573.2444 -12771.0444
58 -8965.2444 -8573.2444
59 -11720.4444 -8965.2444
60 -7697.0556 -11720.4444
61 -4573.5556 -7697.0556
62 2712.1111 -4573.5556
> 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/71ht31229461290.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/8d45q1229461290.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/9pl1v1229461290.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/10pxgj1229461290.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/11qyrf1229461290.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/12kxgq1229461290.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/13ce7c1229461290.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/14o5jy1229461290.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/15cfn71229461290.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/16fzi71229461290.tab")
+ }
>
> system("convert tmp/1xqr91229461290.ps tmp/1xqr91229461290.png")
> system("convert tmp/26w7e1229461290.ps tmp/26w7e1229461290.png")
> system("convert tmp/3bbdx1229461290.ps tmp/3bbdx1229461290.png")
> system("convert tmp/4bo3o1229461290.ps tmp/4bo3o1229461290.png")
> system("convert tmp/5p90s1229461290.ps tmp/5p90s1229461290.png")
> system("convert tmp/6ro2f1229461290.ps tmp/6ro2f1229461290.png")
> system("convert tmp/71ht31229461290.ps tmp/71ht31229461290.png")
> system("convert tmp/8d45q1229461290.ps tmp/8d45q1229461290.png")
> system("convert tmp/9pl1v1229461290.ps tmp/9pl1v1229461290.png")
> system("convert tmp/10pxgj1229461290.ps tmp/10pxgj1229461290.png")
>
>
> proc.time()
user system elapsed
2.568 1.668 5.353