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(9700,0,9081,0,9084,0,9743,0,8587,0,9731,0,9563,0,9998,0,9437,0,10038,0,9918,0,9252,0,9737,0,9035,0,9133,0,9487,0,8700,0,9627,0,8947,0,9283,0,8829,0,9947,0,9628,0,9318,0,9605,0,8640,0,9214,0,9567,0,8547,0,9185,0,9470,0,9123,0,9278,0,10170,0,9434,0,9655,0,9429,0,8739,0,9552,0,9687,0,9019,1,9672,1,9206,1,9069,1,9788,1,10312,1,10105,1,9863,1,9656,1,9295,1,9946,1,9701,1,9049,1,10190,1,9706,1,9765,1,9893,1,9994,1,10433,1,10073,1,10112,1,9266,1,9820,1,10097,1,9115,1,10411,1,9678,1,10408,1,10153,1,10368,1,10581,1,10597,1,10680,1,9738,1,9556,1),dim=c(2,75),dimnames=list(c('Birth','x'),1:75))
> y <- array(NA,dim=c(2,75),dimnames=list(c('Birth','x'),1:75))
> 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
Birth x M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 9700 0 1 0 0 0 0 0 0 0 0 0 0
2 9081 0 0 1 0 0 0 0 0 0 0 0 0
3 9084 0 0 0 1 0 0 0 0 0 0 0 0
4 9743 0 0 0 0 1 0 0 0 0 0 0 0
5 8587 0 0 0 0 0 1 0 0 0 0 0 0
6 9731 0 0 0 0 0 0 1 0 0 0 0 0
7 9563 0 0 0 0 0 0 0 1 0 0 0 0
8 9998 0 0 0 0 0 0 0 0 1 0 0 0
9 9437 0 0 0 0 0 0 0 0 0 1 0 0
10 10038 0 0 0 0 0 0 0 0 0 0 1 0
11 9918 0 0 0 0 0 0 0 0 0 0 0 1
12 9252 0 0 0 0 0 0 0 0 0 0 0 0
13 9737 0 1 0 0 0 0 0 0 0 0 0 0
14 9035 0 0 1 0 0 0 0 0 0 0 0 0
15 9133 0 0 0 1 0 0 0 0 0 0 0 0
16 9487 0 0 0 0 1 0 0 0 0 0 0 0
17 8700 0 0 0 0 0 1 0 0 0 0 0 0
18 9627 0 0 0 0 0 0 1 0 0 0 0 0
19 8947 0 0 0 0 0 0 0 1 0 0 0 0
20 9283 0 0 0 0 0 0 0 0 1 0 0 0
21 8829 0 0 0 0 0 0 0 0 0 1 0 0
22 9947 0 0 0 0 0 0 0 0 0 0 1 0
23 9628 0 0 0 0 0 0 0 0 0 0 0 1
24 9318 0 0 0 0 0 0 0 0 0 0 0 0
25 9605 0 1 0 0 0 0 0 0 0 0 0 0
26 8640 0 0 1 0 0 0 0 0 0 0 0 0
27 9214 0 0 0 1 0 0 0 0 0 0 0 0
28 9567 0 0 0 0 1 0 0 0 0 0 0 0
29 8547 0 0 0 0 0 1 0 0 0 0 0 0
30 9185 0 0 0 0 0 0 1 0 0 0 0 0
31 9470 0 0 0 0 0 0 0 1 0 0 0 0
32 9123 0 0 0 0 0 0 0 0 1 0 0 0
33 9278 0 0 0 0 0 0 0 0 0 1 0 0
34 10170 0 0 0 0 0 0 0 0 0 0 1 0
35 9434 0 0 0 0 0 0 0 0 0 0 0 1
36 9655 0 0 0 0 0 0 0 0 0 0 0 0
37 9429 0 1 0 0 0 0 0 0 0 0 0 0
38 8739 0 0 1 0 0 0 0 0 0 0 0 0
39 9552 0 0 0 1 0 0 0 0 0 0 0 0
40 9687 0 0 0 0 1 0 0 0 0 0 0 0
41 9019 1 0 0 0 0 1 0 0 0 0 0 0
42 9672 1 0 0 0 0 0 1 0 0 0 0 0
43 9206 1 0 0 0 0 0 0 1 0 0 0 0
44 9069 1 0 0 0 0 0 0 0 1 0 0 0
45 9788 1 0 0 0 0 0 0 0 0 1 0 0
46 10312 1 0 0 0 0 0 0 0 0 0 1 0
47 10105 1 0 0 0 0 0 0 0 0 0 0 1
48 9863 1 0 0 0 0 0 0 0 0 0 0 0
49 9656 1 1 0 0 0 0 0 0 0 0 0 0
50 9295 1 0 1 0 0 0 0 0 0 0 0 0
51 9946 1 0 0 1 0 0 0 0 0 0 0 0
52 9701 1 0 0 0 1 0 0 0 0 0 0 0
53 9049 1 0 0 0 0 1 0 0 0 0 0 0
54 10190 1 0 0 0 0 0 1 0 0 0 0 0
55 9706 1 0 0 0 0 0 0 1 0 0 0 0
56 9765 1 0 0 0 0 0 0 0 1 0 0 0
57 9893 1 0 0 0 0 0 0 0 0 1 0 0
58 9994 1 0 0 0 0 0 0 0 0 0 1 0
59 10433 1 0 0 0 0 0 0 0 0 0 0 1
60 10073 1 0 0 0 0 0 0 0 0 0 0 0
61 10112 1 1 0 0 0 0 0 0 0 0 0 0
62 9266 1 0 1 0 0 0 0 0 0 0 0 0
63 9820 1 0 0 1 0 0 0 0 0 0 0 0
64 10097 1 0 0 0 1 0 0 0 0 0 0 0
65 9115 1 0 0 0 0 1 0 0 0 0 0 0
66 10411 1 0 0 0 0 0 1 0 0 0 0 0
67 9678 1 0 0 0 0 0 0 1 0 0 0 0
68 10408 1 0 0 0 0 0 0 0 1 0 0 0
69 10153 1 0 0 0 0 0 0 0 0 1 0 0
70 10368 1 0 0 0 0 0 0 0 0 0 1 0
71 10581 1 0 0 0 0 0 0 0 0 0 0 1
72 10597 1 0 0 0 0 0 0 0 0 0 0 0
73 10680 1 1 0 0 0 0 0 0 0 0 0 0
74 9738 1 0 1 0 0 0 0 0 0 0 0 0
75 9556 1 0 0 1 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) x M1 M2 M3 M4
9548.422 489.156 87.511 -644.632 -285.917 2.193
M5 M6 M7 M8 M9 M10
-956.833 9.667 -364.667 -185.333 -230.000 345.167
M11
223.500
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-783.245 -160.434 -4.589 143.583 634.911
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 9548.422 121.787 78.403 < 2e-16 ***
x 489.156 66.745 7.329 5.79e-10 ***
M1 87.511 159.686 0.548 0.585646
M2 -644.632 159.686 -4.037 0.000151 ***
M3 -285.917 159.686 -1.790 0.078256 .
M4 2.193 166.013 0.013 0.989504
M5 -956.833 165.640 -5.777 2.65e-07 ***
M6 9.667 165.640 0.058 0.953650
M7 -364.667 165.640 -2.202 0.031427 *
M8 -185.333 165.640 -1.119 0.267503
M9 -230.000 165.640 -1.389 0.169937
M10 345.167 165.640 2.084 0.041303 *
M11 223.500 165.640 1.349 0.182144
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 286.9 on 62 degrees of freedom
Multiple R-squared: 0.7275, Adjusted R-squared: 0.6748
F-statistic: 13.79 on 12 and 62 DF, p-value: 2.556e-13
> 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.063038821 0.12607764 0.9369612
[2,] 0.025323760 0.05064752 0.9746762
[3,] 0.009737074 0.01947415 0.9902629
[4,] 0.149017375 0.29803475 0.8509826
[5,] 0.376322611 0.75264522 0.6236774
[6,] 0.510518425 0.97896315 0.4894816
[7,] 0.407513451 0.81502690 0.5924865
[8,] 0.348488696 0.69697739 0.6515113
[9,] 0.274190031 0.54838006 0.7258100
[10,] 0.202970666 0.40594133 0.7970293
[11,] 0.222576848 0.44515370 0.7774232
[12,] 0.163762367 0.32752473 0.8362376
[13,] 0.113492226 0.22698445 0.8865078
[14,] 0.077072972 0.15414594 0.9229270
[15,] 0.113714303 0.22742861 0.8862857
[16,] 0.102195838 0.20439168 0.8978042
[17,] 0.128764138 0.25752828 0.8712359
[18,] 0.095711508 0.19142302 0.9042885
[19,] 0.086717737 0.17343547 0.9132823
[20,] 0.091257071 0.18251414 0.9087429
[21,] 0.082789822 0.16557964 0.9172102
[22,] 0.071502976 0.14300595 0.9284970
[23,] 0.062105909 0.12421182 0.9378941
[24,] 0.060277100 0.12055420 0.9397229
[25,] 0.040047888 0.08009578 0.9599521
[26,] 0.025343772 0.05068754 0.9746562
[27,] 0.030090421 0.06018084 0.9699096
[28,] 0.035198089 0.07039618 0.9648019
[29,] 0.176960483 0.35392097 0.8230395
[30,] 0.203300678 0.40660136 0.7966993
[31,] 0.155090653 0.31018131 0.8449093
[32,] 0.158106322 0.31621264 0.8418937
[33,] 0.177541106 0.35508221 0.8224589
[34,] 0.349503972 0.69900794 0.6504960
[35,] 0.301108125 0.60221625 0.6988919
[36,] 0.305339001 0.61067800 0.6946610
[37,] 0.295185361 0.59037072 0.7048146
[38,] 0.218853315 0.43770663 0.7811467
[39,] 0.193341917 0.38668383 0.8066581
[40,] 0.131213343 0.26242669 0.8687867
[41,] 0.206155424 0.41231085 0.7938446
[42,] 0.163711540 0.32742308 0.8362885
[43,] 0.150323224 0.30064645 0.8496768
[44,] 0.094962593 0.18992519 0.9050374
> postscript(file="/var/www/html/rcomp/tmp/14rmh1291978808.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> plot(x[,1], type='l', main='Actuals and Interpolation', ylab='value of Actuals and Interpolation (dots)', xlab='time or index')
> points(x[,1]-mysum$resid)
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/2f0421291978808.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> plot(mysum$resid, type='b', pch=19, main='Residuals', ylab='value of Residuals', xlab='time or index')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/3f0421291978808.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> hist(mysum$resid, main='Residual Histogram', xlab='values of Residuals')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/4f0421291978808.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> densityplot(~mysum$resid,col='black',main='Residual Density Plot', xlab='values of Residuals')
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/5q9l51291978808.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> qqnorm(mysum$resid, main='Residual Normal Q-Q Plot')
> qqline(mysum$resid)
> grid()
> dev.off()
null device
1
> (myerror <- as.ts(mysum$resid))
Time Series:
Start = 1
End = 75
Frequency = 1
1 2 3 4 5 6
64.066826 177.209683 -178.504602 192.385309 -4.588703 172.911297
7 8 9 10 11 12
379.244631 634.911297 118.577964 144.411297 146.077964 -296.422036
13 14 15 16 17 18
101.066826 131.209683 -129.504602 -63.614691 108.411297 68.911297
19 20 21 22 23 24
-236.755369 -80.088703 -489.422036 53.411297 -143.922036 -230.422036
25 26 27 28 29 30
-30.933174 -263.790317 -48.504602 16.385309 -44.588703 -373.088703
31 32 33 34 35 36
286.244631 -240.088703 -40.422036 276.411297 -337.922036 106.577964
37 38 39 40 41 42
-206.933174 -164.790317 289.495398 136.385309 -61.744631 -375.244631
43 44 45 46 47 48
-466.911297 -783.244631 -19.577964 -70.744631 -156.077964 -174.577964
49 50 51 52 53 54
-469.089102 -97.946244 194.339470 -338.770619 -31.744631 142.755369
55 56 57 58 59 60
33.088703 -87.244631 85.422036 -388.744631 171.922036 35.422036
61 62 63 64 65 66
-13.089102 -126.946244 68.339470 57.229381 34.255369 363.755369
67 68 69 70 71 72
5.088703 555.755369 345.422036 -14.744631 319.922036 559.422036
73 74 75
554.910898 345.053756 -195.660530
> postscript(file="/var/www/html/rcomp/tmp/6q9l51291978808.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> dum <- cbind(lag(myerror,k=1),myerror)
> dum
Time Series:
Start = 0
End = 75
Frequency = 1
lag(myerror, k = 1) myerror
0 64.066826 NA
1 177.209683 64.066826
2 -178.504602 177.209683
3 192.385309 -178.504602
4 -4.588703 192.385309
5 172.911297 -4.588703
6 379.244631 172.911297
7 634.911297 379.244631
8 118.577964 634.911297
9 144.411297 118.577964
10 146.077964 144.411297
11 -296.422036 146.077964
12 101.066826 -296.422036
13 131.209683 101.066826
14 -129.504602 131.209683
15 -63.614691 -129.504602
16 108.411297 -63.614691
17 68.911297 108.411297
18 -236.755369 68.911297
19 -80.088703 -236.755369
20 -489.422036 -80.088703
21 53.411297 -489.422036
22 -143.922036 53.411297
23 -230.422036 -143.922036
24 -30.933174 -230.422036
25 -263.790317 -30.933174
26 -48.504602 -263.790317
27 16.385309 -48.504602
28 -44.588703 16.385309
29 -373.088703 -44.588703
30 286.244631 -373.088703
31 -240.088703 286.244631
32 -40.422036 -240.088703
33 276.411297 -40.422036
34 -337.922036 276.411297
35 106.577964 -337.922036
36 -206.933174 106.577964
37 -164.790317 -206.933174
38 289.495398 -164.790317
39 136.385309 289.495398
40 -61.744631 136.385309
41 -375.244631 -61.744631
42 -466.911297 -375.244631
43 -783.244631 -466.911297
44 -19.577964 -783.244631
45 -70.744631 -19.577964
46 -156.077964 -70.744631
47 -174.577964 -156.077964
48 -469.089102 -174.577964
49 -97.946244 -469.089102
50 194.339470 -97.946244
51 -338.770619 194.339470
52 -31.744631 -338.770619
53 142.755369 -31.744631
54 33.088703 142.755369
55 -87.244631 33.088703
56 85.422036 -87.244631
57 -388.744631 85.422036
58 171.922036 -388.744631
59 35.422036 171.922036
60 -13.089102 35.422036
61 -126.946244 -13.089102
62 68.339470 -126.946244
63 57.229381 68.339470
64 34.255369 57.229381
65 363.755369 34.255369
66 5.088703 363.755369
67 555.755369 5.088703
68 345.422036 555.755369
69 -14.744631 345.422036
70 319.922036 -14.744631
71 559.422036 319.922036
72 554.910898 559.422036
73 345.053756 554.910898
74 -195.660530 345.053756
75 NA -195.660530
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 177.209683 64.066826
[2,] -178.504602 177.209683
[3,] 192.385309 -178.504602
[4,] -4.588703 192.385309
[5,] 172.911297 -4.588703
[6,] 379.244631 172.911297
[7,] 634.911297 379.244631
[8,] 118.577964 634.911297
[9,] 144.411297 118.577964
[10,] 146.077964 144.411297
[11,] -296.422036 146.077964
[12,] 101.066826 -296.422036
[13,] 131.209683 101.066826
[14,] -129.504602 131.209683
[15,] -63.614691 -129.504602
[16,] 108.411297 -63.614691
[17,] 68.911297 108.411297
[18,] -236.755369 68.911297
[19,] -80.088703 -236.755369
[20,] -489.422036 -80.088703
[21,] 53.411297 -489.422036
[22,] -143.922036 53.411297
[23,] -230.422036 -143.922036
[24,] -30.933174 -230.422036
[25,] -263.790317 -30.933174
[26,] -48.504602 -263.790317
[27,] 16.385309 -48.504602
[28,] -44.588703 16.385309
[29,] -373.088703 -44.588703
[30,] 286.244631 -373.088703
[31,] -240.088703 286.244631
[32,] -40.422036 -240.088703
[33,] 276.411297 -40.422036
[34,] -337.922036 276.411297
[35,] 106.577964 -337.922036
[36,] -206.933174 106.577964
[37,] -164.790317 -206.933174
[38,] 289.495398 -164.790317
[39,] 136.385309 289.495398
[40,] -61.744631 136.385309
[41,] -375.244631 -61.744631
[42,] -466.911297 -375.244631
[43,] -783.244631 -466.911297
[44,] -19.577964 -783.244631
[45,] -70.744631 -19.577964
[46,] -156.077964 -70.744631
[47,] -174.577964 -156.077964
[48,] -469.089102 -174.577964
[49,] -97.946244 -469.089102
[50,] 194.339470 -97.946244
[51,] -338.770619 194.339470
[52,] -31.744631 -338.770619
[53,] 142.755369 -31.744631
[54,] 33.088703 142.755369
[55,] -87.244631 33.088703
[56,] 85.422036 -87.244631
[57,] -388.744631 85.422036
[58,] 171.922036 -388.744631
[59,] 35.422036 171.922036
[60,] -13.089102 35.422036
[61,] -126.946244 -13.089102
[62,] 68.339470 -126.946244
[63,] 57.229381 68.339470
[64,] 34.255369 57.229381
[65,] 363.755369 34.255369
[66,] 5.088703 363.755369
[67,] 555.755369 5.088703
[68,] 345.422036 555.755369
[69,] -14.744631 345.422036
[70,] 319.922036 -14.744631
[71,] 559.422036 319.922036
[72,] 554.910898 559.422036
[73,] 345.053756 554.910898
[74,] -195.660530 345.053756
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 177.209683 64.066826
2 -178.504602 177.209683
3 192.385309 -178.504602
4 -4.588703 192.385309
5 172.911297 -4.588703
6 379.244631 172.911297
7 634.911297 379.244631
8 118.577964 634.911297
9 144.411297 118.577964
10 146.077964 144.411297
11 -296.422036 146.077964
12 101.066826 -296.422036
13 131.209683 101.066826
14 -129.504602 131.209683
15 -63.614691 -129.504602
16 108.411297 -63.614691
17 68.911297 108.411297
18 -236.755369 68.911297
19 -80.088703 -236.755369
20 -489.422036 -80.088703
21 53.411297 -489.422036
22 -143.922036 53.411297
23 -230.422036 -143.922036
24 -30.933174 -230.422036
25 -263.790317 -30.933174
26 -48.504602 -263.790317
27 16.385309 -48.504602
28 -44.588703 16.385309
29 -373.088703 -44.588703
30 286.244631 -373.088703
31 -240.088703 286.244631
32 -40.422036 -240.088703
33 276.411297 -40.422036
34 -337.922036 276.411297
35 106.577964 -337.922036
36 -206.933174 106.577964
37 -164.790317 -206.933174
38 289.495398 -164.790317
39 136.385309 289.495398
40 -61.744631 136.385309
41 -375.244631 -61.744631
42 -466.911297 -375.244631
43 -783.244631 -466.911297
44 -19.577964 -783.244631
45 -70.744631 -19.577964
46 -156.077964 -70.744631
47 -174.577964 -156.077964
48 -469.089102 -174.577964
49 -97.946244 -469.089102
50 194.339470 -97.946244
51 -338.770619 194.339470
52 -31.744631 -338.770619
53 142.755369 -31.744631
54 33.088703 142.755369
55 -87.244631 33.088703
56 85.422036 -87.244631
57 -388.744631 85.422036
58 171.922036 -388.744631
59 35.422036 171.922036
60 -13.089102 35.422036
61 -126.946244 -13.089102
62 68.339470 -126.946244
63 57.229381 68.339470
64 34.255369 57.229381
65 363.755369 34.255369
66 5.088703 363.755369
67 555.755369 5.088703
68 345.422036 555.755369
69 -14.744631 345.422036
70 319.922036 -14.744631
71 559.422036 319.922036
72 554.910898 559.422036
73 345.053756 554.910898
74 -195.660530 345.053756
> 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/7i1k81291978808.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> acf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Autocorrelation Function')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/8bs1t1291978808.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> pacf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Partial Autocorrelation Function')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/9bs1t1291978808.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> opar <- par(mfrow = c(2,2), oma = c(0, 0, 1.1, 0))
> plot(mylm, las = 1, sub='Residual Diagnostics')
> par(opar)
> dev.off()
null device
1
> if (n > n25) {
+ postscript(file="/var/www/html/rcomp/tmp/10bs1t1291978808.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
+ plot(kp3:nmkm3,gqarr[,2], main='Goldfeld-Quandt test',ylab='2-sided p-value',xlab='breakpoint')
+ grid()
+ dev.off()
+ }
null device
1
>
> #Note: the /var/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/11p2zk1291978808.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/12bkyp1291978808.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/13h4yt1291978809.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/14k4ez1291978809.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/15vwek1291978809.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/16r5ts1291978809.tab")
+ }
>
> try(system("convert tmp/14rmh1291978808.ps tmp/14rmh1291978808.png",intern=TRUE))
character(0)
> try(system("convert tmp/2f0421291978808.ps tmp/2f0421291978808.png",intern=TRUE))
character(0)
> try(system("convert tmp/3f0421291978808.ps tmp/3f0421291978808.png",intern=TRUE))
character(0)
> try(system("convert tmp/4f0421291978808.ps tmp/4f0421291978808.png",intern=TRUE))
character(0)
> try(system("convert tmp/5q9l51291978808.ps tmp/5q9l51291978808.png",intern=TRUE))
character(0)
> try(system("convert tmp/6q9l51291978808.ps tmp/6q9l51291978808.png",intern=TRUE))
character(0)
> try(system("convert tmp/7i1k81291978808.ps tmp/7i1k81291978808.png",intern=TRUE))
character(0)
> try(system("convert tmp/8bs1t1291978808.ps tmp/8bs1t1291978808.png",intern=TRUE))
character(0)
> try(system("convert tmp/9bs1t1291978808.ps tmp/9bs1t1291978808.png",intern=TRUE))
character(0)
> try(system("convert tmp/10bs1t1291978808.ps tmp/10bs1t1291978808.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
2.660 1.648 8.778