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.
Natural language support but running in an English locale
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(9081,0,9700,9084,0,9081,9743,0,9084,8587,0,9743,9731,0,8587,9563,0,9731,9998,0,9563,9437,0,9998,10038,0,9437,9918,0,10038,9252,0,9918,9737,0,9252,9035,0,9737,9133,0,9035,9487,0,9133,8700,0,9487,9627,0,8700,8947,0,9627,9283,0,8947,8829,0,9283,9947,0,8829,9628,0,9947,9318,0,9628,9605,0,9318,8640,0,9605,9214,0,8640,9567,0,9214,8547,0,9567,9185,0,8547,9470,0,9185,9123,0,9470,9278,0,9123,10170,0,9278,9434,0,10170,9655,0,9434,9429,0,9655,8739,0,9429,9552,0,8739,9687,1,9552,9019,1,9687,9672,1,9019,9206,1,9672,9069,1,9206,9788,1,9069,10312,1,9788,10105,1,10312,9863,1,10105,9656,1,9863,9295,1,9656,9946,1,9295,9701,1,9946,9049,1,9701,10190,1,9049,9706,1,10190,9765,1,9706,9893,1,9765,9994,1,9893,10433,1,9994,10073,1,10433,10112,1,10073,9266,1,10112,9820,1,9266,10097,1,9820,9115,1,10097,10411,1,9115,9678,1,10411,10408,1,9678,10153,1,10408,10368,1,10153,10581,1,10368,10597,1,10581,10680,1,10597,9738,1,10680,9556,1,9738),dim=c(3,74),dimnames=list(c('geboortes','x','lag'),1:74))
> y <- array(NA,dim=c(3,74),dimnames=list(c('geboortes','x','lag'),1:74))
> 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
geboortes x lag M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 9081 0 9700 1 0 0 0 0 0 0 0 0 0 0 1
2 9084 0 9081 0 1 0 0 0 0 0 0 0 0 0 2
3 9743 0 9084 0 0 1 0 0 0 0 0 0 0 0 3
4 8587 0 9743 0 0 0 1 0 0 0 0 0 0 0 4
5 9731 0 8587 0 0 0 0 1 0 0 0 0 0 0 5
6 9563 0 9731 0 0 0 0 0 1 0 0 0 0 0 6
7 9998 0 9563 0 0 0 0 0 0 1 0 0 0 0 7
8 9437 0 9998 0 0 0 0 0 0 0 1 0 0 0 8
9 10038 0 9437 0 0 0 0 0 0 0 0 1 0 0 9
10 9918 0 10038 0 0 0 0 0 0 0 0 0 1 0 10
11 9252 0 9918 0 0 0 0 0 0 0 0 0 0 1 11
12 9737 0 9252 0 0 0 0 0 0 0 0 0 0 0 12
13 9035 0 9737 1 0 0 0 0 0 0 0 0 0 0 13
14 9133 0 9035 0 1 0 0 0 0 0 0 0 0 0 14
15 9487 0 9133 0 0 1 0 0 0 0 0 0 0 0 15
16 8700 0 9487 0 0 0 1 0 0 0 0 0 0 0 16
17 9627 0 8700 0 0 0 0 1 0 0 0 0 0 0 17
18 8947 0 9627 0 0 0 0 0 1 0 0 0 0 0 18
19 9283 0 8947 0 0 0 0 0 0 1 0 0 0 0 19
20 8829 0 9283 0 0 0 0 0 0 0 1 0 0 0 20
21 9947 0 8829 0 0 0 0 0 0 0 0 1 0 0 21
22 9628 0 9947 0 0 0 0 0 0 0 0 0 1 0 22
23 9318 0 9628 0 0 0 0 0 0 0 0 0 0 1 23
24 9605 0 9318 0 0 0 0 0 0 0 0 0 0 0 24
25 8640 0 9605 1 0 0 0 0 0 0 0 0 0 0 25
26 9214 0 8640 0 1 0 0 0 0 0 0 0 0 0 26
27 9567 0 9214 0 0 1 0 0 0 0 0 0 0 0 27
28 8547 0 9567 0 0 0 1 0 0 0 0 0 0 0 28
29 9185 0 8547 0 0 0 0 1 0 0 0 0 0 0 29
30 9470 0 9185 0 0 0 0 0 1 0 0 0 0 0 30
31 9123 0 9470 0 0 0 0 0 0 1 0 0 0 0 31
32 9278 0 9123 0 0 0 0 0 0 0 1 0 0 0 32
33 10170 0 9278 0 0 0 0 0 0 0 0 1 0 0 33
34 9434 0 10170 0 0 0 0 0 0 0 0 0 1 0 34
35 9655 0 9434 0 0 0 0 0 0 0 0 0 0 1 35
36 9429 0 9655 0 0 0 0 0 0 0 0 0 0 0 36
37 8739 0 9429 1 0 0 0 0 0 0 0 0 0 0 37
38 9552 0 8739 0 1 0 0 0 0 0 0 0 0 0 38
39 9687 1 9552 0 0 1 0 0 0 0 0 0 0 0 39
40 9019 1 9687 0 0 0 1 0 0 0 0 0 0 0 40
41 9672 1 9019 0 0 0 0 1 0 0 0 0 0 0 41
42 9206 1 9672 0 0 0 0 0 1 0 0 0 0 0 42
43 9069 1 9206 0 0 0 0 0 0 1 0 0 0 0 43
44 9788 1 9069 0 0 0 0 0 0 0 1 0 0 0 44
45 10312 1 9788 0 0 0 0 0 0 0 0 1 0 0 45
46 10105 1 10312 0 0 0 0 0 0 0 0 0 1 0 46
47 9863 1 10105 0 0 0 0 0 0 0 0 0 0 1 47
48 9656 1 9863 0 0 0 0 0 0 0 0 0 0 0 48
49 9295 1 9656 1 0 0 0 0 0 0 0 0 0 0 49
50 9946 1 9295 0 1 0 0 0 0 0 0 0 0 0 50
51 9701 1 9946 0 0 1 0 0 0 0 0 0 0 0 51
52 9049 1 9701 0 0 0 1 0 0 0 0 0 0 0 52
53 10190 1 9049 0 0 0 0 1 0 0 0 0 0 0 53
54 9706 1 10190 0 0 0 0 0 1 0 0 0 0 0 54
55 9765 1 9706 0 0 0 0 0 0 1 0 0 0 0 55
56 9893 1 9765 0 0 0 0 0 0 0 1 0 0 0 56
57 9994 1 9893 0 0 0 0 0 0 0 0 1 0 0 57
58 10433 1 9994 0 0 0 0 0 0 0 0 0 1 0 58
59 10073 1 10433 0 0 0 0 0 0 0 0 0 0 1 59
60 10112 1 10073 0 0 0 0 0 0 0 0 0 0 0 60
61 9266 1 10112 1 0 0 0 0 0 0 0 0 0 0 61
62 9820 1 9266 0 1 0 0 0 0 0 0 0 0 0 62
63 10097 1 9820 0 0 1 0 0 0 0 0 0 0 0 63
64 9115 1 10097 0 0 0 1 0 0 0 0 0 0 0 64
65 10411 1 9115 0 0 0 0 1 0 0 0 0 0 0 65
66 9678 1 10411 0 0 0 0 0 1 0 0 0 0 0 66
67 10408 1 9678 0 0 0 0 0 0 1 0 0 0 0 67
68 10153 1 10408 0 0 0 0 0 0 0 1 0 0 0 68
69 10368 1 10153 0 0 0 0 0 0 0 0 1 0 0 69
70 10581 1 10368 0 0 0 0 0 0 0 0 0 1 0 70
71 10597 1 10581 0 0 0 0 0 0 0 0 0 0 1 71
72 10680 1 10597 0 0 0 0 0 0 0 0 0 0 0 72
73 9738 1 10680 1 0 0 0 0 0 0 0 0 0 0 73
74 9556 1 9738 0 1 0 0 0 0 0 0 0 0 0 74
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) x lag M1 M2 M3
7081.7941 199.6758 0.2561 -734.1016 -192.2160 -31.7261
M4 M5 M6 M7 M8 M9
-978.9495 207.9418 -418.1729 -147.2886 -242.1755 340.1281
M10 M11 t
66.8845 -129.7621 4.3004
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-607.377 -159.252 2.071 159.764 584.700
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 7081.7941 1203.0951 5.886 2.00e-07 ***
x 199.6758 138.3700 1.443 0.1543
lag 0.2561 0.1269 2.018 0.0481 *
M1 -734.1016 155.9461 -4.707 1.56e-05 ***
M2 -192.2160 175.0807 -1.098 0.2767
M3 -31.7261 167.4559 -0.189 0.8504
M4 -978.9495 163.0476 -6.004 1.27e-07 ***
M5 207.9418 199.9725 1.040 0.3027
M6 -418.1729 162.2269 -2.578 0.0125 *
M7 -147.2886 167.3572 -0.880 0.3824
M8 -242.1755 162.8483 -1.487 0.1423
M9 340.1281 163.5689 2.079 0.0419 *
M10 66.8845 167.3761 0.400 0.6909
M11 -129.7621 163.6667 -0.793 0.4310
t 4.3004 3.2199 1.336 0.1868
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 278.9 on 59 degrees of freedom
Multiple R-squared: 0.7548, Adjusted R-squared: 0.6966
F-statistic: 12.97 on 14 and 59 DF, p-value: 3.444e-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.6010485 0.7979031 0.3989515
[2,] 0.5767998 0.8464004 0.4232002
[3,] 0.4408215 0.8816430 0.5591785
[4,] 0.4889042 0.9778084 0.5110958
[5,] 0.3795509 0.7591018 0.6204491
[6,] 0.3710491 0.7420982 0.6289509
[7,] 0.2802361 0.5604722 0.7197639
[8,] 0.2151753 0.4303506 0.7848247
[9,] 0.3222142 0.6444284 0.6777858
[10,] 0.2580855 0.5161710 0.7419145
[11,] 0.1892858 0.3785716 0.8107142
[12,] 0.2089239 0.4178478 0.7910761
[13,] 0.3593077 0.7186155 0.6406923
[14,] 0.3526087 0.7052173 0.6473913
[15,] 0.3555199 0.7110398 0.6444801
[16,] 0.4493296 0.8986592 0.5506704
[17,] 0.4218910 0.8437820 0.5781090
[18,] 0.4878443 0.9756887 0.5121557
[19,] 0.4244597 0.8489193 0.5755403
[20,] 0.3932366 0.7864733 0.6067634
[21,] 0.4787794 0.9575589 0.5212206
[22,] 0.4034393 0.8068786 0.5965607
[23,] 0.3852674 0.7705347 0.6147326
[24,] 0.3309876 0.6619753 0.6690124
[25,] 0.2950761 0.5901522 0.7049239
[26,] 0.5507747 0.8984505 0.4492253
[27,] 0.5511121 0.8977759 0.4488879
[28,] 0.6150218 0.7699565 0.3849782
[29,] 0.5478055 0.9043890 0.4521945
[30,] 0.4837766 0.9675532 0.5162234
[31,] 0.5313080 0.9373840 0.4686920
[32,] 0.4481055 0.8962110 0.5518945
[33,] 0.7771815 0.4456370 0.2228185
[34,] 0.6814095 0.6371810 0.3185905
[35,] 0.6084117 0.7831767 0.3915883
[36,] 0.5496218 0.9007564 0.4503782
[37,] 0.5812686 0.8374629 0.4187314
[38,] 0.4580159 0.9160318 0.5419841
[39,] 0.3137093 0.6274186 0.6862907
> postscript(file="/var/www/html/freestat/rcomp/tmp/18zys1291970973.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/freestat/rcomp/tmp/2i8fv1291970973.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/freestat/rcomp/tmp/3i8fv1291970973.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/freestat/rcomp/tmp/4i8fv1291970973.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/freestat/rcomp/tmp/5t0fy1291970973.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 = 74
Frequency = 1
1 2 3 4 5 6
245.2357663 -139.4497553 353.9917240 -27.8280435 220.9843610 381.8672391
7 8 9 10 11 12
584.7004187 2.9013591 160.9464395 155.9982511 -286.9285137 234.5442177
13 14 15 16 17 18
138.1568809 -130.2757517 33.8401317 99.1183298 36.4449990 -259.1073410
19 20 21 22 23 24
-24.1720029 -473.6212310 174.0255467 -162.3050946 -198.2761377 34.0396241
25 26 27 28 29 30
-274.6480498 0.2628091 41.4946546 -125.9710884 -417.9826947 325.4659883
31 32 33 34 35 36
-369.6955152 -35.2565125 230.4503932 -465.0109360 136.7945836 -279.8569323
37 38 39 40 41 42
-182.1863888 261.3082717 -176.3337995 64.0212983 -303.1230419 -314.5152424
43 44 45 46 47 48
-607.3765096 237.2901235 -9.4201923 -81.6518451 -78.3014854 -357.3977290
49 50 51 52 53 54
64.4076954 261.6589767 -314.8257133 38.8317677 155.5904851 1.2415398
55 56 57 58 59 60
-91.0106672 112.4684209 -405.9110831 276.1701802 -3.8935116 -6.7748046
61 62 63 64 65 66
-132.9598705 91.4799789 61.8330025 -48.1722639 308.0858915 -134.9521838
67 68 69 70 71 72
507.5542762 156.2178401 -150.0911039 276.7994444 430.6050648 375.4456241
73 74
141.9939665 -344.9845294
> postscript(file="/var/www/html/freestat/rcomp/tmp/6t0fy1291970973.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 = 74
Frequency = 1
lag(myerror, k = 1) myerror
0 245.2357663 NA
1 -139.4497553 245.2357663
2 353.9917240 -139.4497553
3 -27.8280435 353.9917240
4 220.9843610 -27.8280435
5 381.8672391 220.9843610
6 584.7004187 381.8672391
7 2.9013591 584.7004187
8 160.9464395 2.9013591
9 155.9982511 160.9464395
10 -286.9285137 155.9982511
11 234.5442177 -286.9285137
12 138.1568809 234.5442177
13 -130.2757517 138.1568809
14 33.8401317 -130.2757517
15 99.1183298 33.8401317
16 36.4449990 99.1183298
17 -259.1073410 36.4449990
18 -24.1720029 -259.1073410
19 -473.6212310 -24.1720029
20 174.0255467 -473.6212310
21 -162.3050946 174.0255467
22 -198.2761377 -162.3050946
23 34.0396241 -198.2761377
24 -274.6480498 34.0396241
25 0.2628091 -274.6480498
26 41.4946546 0.2628091
27 -125.9710884 41.4946546
28 -417.9826947 -125.9710884
29 325.4659883 -417.9826947
30 -369.6955152 325.4659883
31 -35.2565125 -369.6955152
32 230.4503932 -35.2565125
33 -465.0109360 230.4503932
34 136.7945836 -465.0109360
35 -279.8569323 136.7945836
36 -182.1863888 -279.8569323
37 261.3082717 -182.1863888
38 -176.3337995 261.3082717
39 64.0212983 -176.3337995
40 -303.1230419 64.0212983
41 -314.5152424 -303.1230419
42 -607.3765096 -314.5152424
43 237.2901235 -607.3765096
44 -9.4201923 237.2901235
45 -81.6518451 -9.4201923
46 -78.3014854 -81.6518451
47 -357.3977290 -78.3014854
48 64.4076954 -357.3977290
49 261.6589767 64.4076954
50 -314.8257133 261.6589767
51 38.8317677 -314.8257133
52 155.5904851 38.8317677
53 1.2415398 155.5904851
54 -91.0106672 1.2415398
55 112.4684209 -91.0106672
56 -405.9110831 112.4684209
57 276.1701802 -405.9110831
58 -3.8935116 276.1701802
59 -6.7748046 -3.8935116
60 -132.9598705 -6.7748046
61 91.4799789 -132.9598705
62 61.8330025 91.4799789
63 -48.1722639 61.8330025
64 308.0858915 -48.1722639
65 -134.9521838 308.0858915
66 507.5542762 -134.9521838
67 156.2178401 507.5542762
68 -150.0911039 156.2178401
69 276.7994444 -150.0911039
70 430.6050648 276.7994444
71 375.4456241 430.6050648
72 141.9939665 375.4456241
73 -344.9845294 141.9939665
74 NA -344.9845294
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -139.4497553 245.2357663
[2,] 353.9917240 -139.4497553
[3,] -27.8280435 353.9917240
[4,] 220.9843610 -27.8280435
[5,] 381.8672391 220.9843610
[6,] 584.7004187 381.8672391
[7,] 2.9013591 584.7004187
[8,] 160.9464395 2.9013591
[9,] 155.9982511 160.9464395
[10,] -286.9285137 155.9982511
[11,] 234.5442177 -286.9285137
[12,] 138.1568809 234.5442177
[13,] -130.2757517 138.1568809
[14,] 33.8401317 -130.2757517
[15,] 99.1183298 33.8401317
[16,] 36.4449990 99.1183298
[17,] -259.1073410 36.4449990
[18,] -24.1720029 -259.1073410
[19,] -473.6212310 -24.1720029
[20,] 174.0255467 -473.6212310
[21,] -162.3050946 174.0255467
[22,] -198.2761377 -162.3050946
[23,] 34.0396241 -198.2761377
[24,] -274.6480498 34.0396241
[25,] 0.2628091 -274.6480498
[26,] 41.4946546 0.2628091
[27,] -125.9710884 41.4946546
[28,] -417.9826947 -125.9710884
[29,] 325.4659883 -417.9826947
[30,] -369.6955152 325.4659883
[31,] -35.2565125 -369.6955152
[32,] 230.4503932 -35.2565125
[33,] -465.0109360 230.4503932
[34,] 136.7945836 -465.0109360
[35,] -279.8569323 136.7945836
[36,] -182.1863888 -279.8569323
[37,] 261.3082717 -182.1863888
[38,] -176.3337995 261.3082717
[39,] 64.0212983 -176.3337995
[40,] -303.1230419 64.0212983
[41,] -314.5152424 -303.1230419
[42,] -607.3765096 -314.5152424
[43,] 237.2901235 -607.3765096
[44,] -9.4201923 237.2901235
[45,] -81.6518451 -9.4201923
[46,] -78.3014854 -81.6518451
[47,] -357.3977290 -78.3014854
[48,] 64.4076954 -357.3977290
[49,] 261.6589767 64.4076954
[50,] -314.8257133 261.6589767
[51,] 38.8317677 -314.8257133
[52,] 155.5904851 38.8317677
[53,] 1.2415398 155.5904851
[54,] -91.0106672 1.2415398
[55,] 112.4684209 -91.0106672
[56,] -405.9110831 112.4684209
[57,] 276.1701802 -405.9110831
[58,] -3.8935116 276.1701802
[59,] -6.7748046 -3.8935116
[60,] -132.9598705 -6.7748046
[61,] 91.4799789 -132.9598705
[62,] 61.8330025 91.4799789
[63,] -48.1722639 61.8330025
[64,] 308.0858915 -48.1722639
[65,] -134.9521838 308.0858915
[66,] 507.5542762 -134.9521838
[67,] 156.2178401 507.5542762
[68,] -150.0911039 156.2178401
[69,] 276.7994444 -150.0911039
[70,] 430.6050648 276.7994444
[71,] 375.4456241 430.6050648
[72,] 141.9939665 375.4456241
[73,] -344.9845294 141.9939665
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -139.4497553 245.2357663
2 353.9917240 -139.4497553
3 -27.8280435 353.9917240
4 220.9843610 -27.8280435
5 381.8672391 220.9843610
6 584.7004187 381.8672391
7 2.9013591 584.7004187
8 160.9464395 2.9013591
9 155.9982511 160.9464395
10 -286.9285137 155.9982511
11 234.5442177 -286.9285137
12 138.1568809 234.5442177
13 -130.2757517 138.1568809
14 33.8401317 -130.2757517
15 99.1183298 33.8401317
16 36.4449990 99.1183298
17 -259.1073410 36.4449990
18 -24.1720029 -259.1073410
19 -473.6212310 -24.1720029
20 174.0255467 -473.6212310
21 -162.3050946 174.0255467
22 -198.2761377 -162.3050946
23 34.0396241 -198.2761377
24 -274.6480498 34.0396241
25 0.2628091 -274.6480498
26 41.4946546 0.2628091
27 -125.9710884 41.4946546
28 -417.9826947 -125.9710884
29 325.4659883 -417.9826947
30 -369.6955152 325.4659883
31 -35.2565125 -369.6955152
32 230.4503932 -35.2565125
33 -465.0109360 230.4503932
34 136.7945836 -465.0109360
35 -279.8569323 136.7945836
36 -182.1863888 -279.8569323
37 261.3082717 -182.1863888
38 -176.3337995 261.3082717
39 64.0212983 -176.3337995
40 -303.1230419 64.0212983
41 -314.5152424 -303.1230419
42 -607.3765096 -314.5152424
43 237.2901235 -607.3765096
44 -9.4201923 237.2901235
45 -81.6518451 -9.4201923
46 -78.3014854 -81.6518451
47 -357.3977290 -78.3014854
48 64.4076954 -357.3977290
49 261.6589767 64.4076954
50 -314.8257133 261.6589767
51 38.8317677 -314.8257133
52 155.5904851 38.8317677
53 1.2415398 155.5904851
54 -91.0106672 1.2415398
55 112.4684209 -91.0106672
56 -405.9110831 112.4684209
57 276.1701802 -405.9110831
58 -3.8935116 276.1701802
59 -6.7748046 -3.8935116
60 -132.9598705 -6.7748046
61 91.4799789 -132.9598705
62 61.8330025 91.4799789
63 -48.1722639 61.8330025
64 308.0858915 -48.1722639
65 -134.9521838 308.0858915
66 507.5542762 -134.9521838
67 156.2178401 507.5542762
68 -150.0911039 156.2178401
69 276.7994444 -150.0911039
70 430.6050648 276.7994444
71 375.4456241 430.6050648
72 141.9939665 375.4456241
73 -344.9845294 141.9939665
> 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/freestat/rcomp/tmp/74rej1291970973.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/freestat/rcomp/tmp/84rej1291970973.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/freestat/rcomp/tmp/9w0d41291970973.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/freestat/rcomp/tmp/10w0d41291970973.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/freestat/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/html/freestat/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/freestat/rcomp/tmp/11ijca1291970973.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/freestat/rcomp/tmp/12ljsy1291970973.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/freestat/rcomp/tmp/13hbq61291970973.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/freestat/rcomp/tmp/14vm971291970974.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/freestat/rcomp/tmp/15hm8d1291970974.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/freestat/rcomp/tmp/16izca1291970974.tab")
+ }
>
> try(system("convert tmp/18zys1291970973.ps tmp/18zys1291970973.png",intern=TRUE))
character(0)
> try(system("convert tmp/2i8fv1291970973.ps tmp/2i8fv1291970973.png",intern=TRUE))
character(0)
> try(system("convert tmp/3i8fv1291970973.ps tmp/3i8fv1291970973.png",intern=TRUE))
character(0)
> try(system("convert tmp/4i8fv1291970973.ps tmp/4i8fv1291970973.png",intern=TRUE))
character(0)
> try(system("convert tmp/5t0fy1291970973.ps tmp/5t0fy1291970973.png",intern=TRUE))
character(0)
> try(system("convert tmp/6t0fy1291970973.ps tmp/6t0fy1291970973.png",intern=TRUE))
character(0)
> try(system("convert tmp/74rej1291970973.ps tmp/74rej1291970973.png",intern=TRUE))
character(0)
> try(system("convert tmp/84rej1291970973.ps tmp/84rej1291970973.png",intern=TRUE))
character(0)
> try(system("convert tmp/9w0d41291970973.ps tmp/9w0d41291970973.png",intern=TRUE))
character(0)
> try(system("convert tmp/10w0d41291970973.ps tmp/10w0d41291970973.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
4.089 2.588 4.464