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(16198.9
+ ,16896.2
+ ,0
+ ,16554.2
+ ,16698
+ ,0
+ ,19554.2
+ ,19691.6
+ ,0
+ ,15903.8
+ ,15930.7
+ ,0
+ ,18003.8
+ ,17444.6
+ ,0
+ ,18329.6
+ ,17699.4
+ ,0
+ ,16260.7
+ ,15189.8
+ ,0
+ ,14851.9
+ ,15672.7
+ ,0
+ ,18174.1
+ ,17180.8
+ ,0
+ ,18406.6
+ ,17664.9
+ ,0
+ ,18466.5
+ ,17862.9
+ ,0
+ ,16016.5
+ ,16162.3
+ ,0
+ ,17428.5
+ ,17463.6
+ ,0
+ ,17167.2
+ ,16772.1
+ ,0
+ ,19630
+ ,19106.9
+ ,0
+ ,17183.6
+ ,16721.3
+ ,0
+ ,18344.7
+ ,18161.3
+ ,0
+ ,19301.4
+ ,18509.9
+ ,0
+ ,18147.5
+ ,17802.7
+ ,0
+ ,16192.9
+ ,16409.9
+ ,0
+ ,18374.4
+ ,17967.7
+ ,0
+ ,20515.2
+ ,20286.6
+ ,0
+ ,18957.2
+ ,19537.3
+ ,0
+ ,16471.5
+ ,18021.9
+ ,0
+ ,18746.8
+ ,20194.3
+ ,0
+ ,19009.5
+ ,19049.6
+ ,0
+ ,19211.2
+ ,20244.7
+ ,0
+ ,20547.7
+ ,21473.3
+ ,0
+ ,19325.8
+ ,19673.6
+ ,0
+ ,20605.5
+ ,21053.2
+ ,0
+ ,20056.9
+ ,20159.5
+ ,0
+ ,16141.4
+ ,18203.6
+ ,0
+ ,20359.8
+ ,21289.5
+ ,0
+ ,19711.6
+ ,20432.3
+ ,1
+ ,15638.6
+ ,17180.4
+ ,1
+ ,14384.5
+ ,15816.8
+ ,1
+ ,13855.6
+ ,15071.8
+ ,1
+ ,14308.3
+ ,14521.1
+ ,1
+ ,15290.6
+ ,15668.8
+ ,1
+ ,14423.8
+ ,14346.9
+ ,1
+ ,13779.7
+ ,13881
+ ,1
+ ,15686.3
+ ,15465.9
+ ,1
+ ,14733.8
+ ,14238.2
+ ,1
+ ,12522.5
+ ,13557.7
+ ,1
+ ,16189.4
+ ,16127.6
+ ,1
+ ,16059.1
+ ,16793.9
+ ,1
+ ,16007.1
+ ,16014
+ ,1
+ ,15806.8
+ ,16867.9
+ ,1
+ ,15160
+ ,16014.6
+ ,0
+ ,15692.1
+ ,15878.6
+ ,0
+ ,18908.9
+ ,18664.9
+ ,0
+ ,16969.9
+ ,17962.5
+ ,0
+ ,16997.5
+ ,17332.7
+ ,0
+ ,19858.9
+ ,19542.1
+ ,0
+ ,17681.2
+ ,17203.6
+ ,0)
+ ,dim=c(3
+ ,55)
+ ,dimnames=list(c('uitvoer'
+ ,'invoer'
+ ,'crisis')
+ ,1:55))
> y <- array(NA,dim=c(3,55),dimnames=list(c('uitvoer','invoer','crisis'),1:55))
> 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
uitvoer invoer crisis M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 16198.9 16896.2 0 1 0 0 0 0 0 0 0 0 0 0 1
2 16554.2 16698.0 0 0 1 0 0 0 0 0 0 0 0 0 2
3 19554.2 19691.6 0 0 0 1 0 0 0 0 0 0 0 0 3
4 15903.8 15930.7 0 0 0 0 1 0 0 0 0 0 0 0 4
5 18003.8 17444.6 0 0 0 0 0 1 0 0 0 0 0 0 5
6 18329.6 17699.4 0 0 0 0 0 0 1 0 0 0 0 0 6
7 16260.7 15189.8 0 0 0 0 0 0 0 1 0 0 0 0 7
8 14851.9 15672.7 0 0 0 0 0 0 0 0 1 0 0 0 8
9 18174.1 17180.8 0 0 0 0 0 0 0 0 0 1 0 0 9
10 18406.6 17664.9 0 0 0 0 0 0 0 0 0 0 1 0 10
11 18466.5 17862.9 0 0 0 0 0 0 0 0 0 0 0 1 11
12 16016.5 16162.3 0 0 0 0 0 0 0 0 0 0 0 0 12
13 17428.5 17463.6 0 1 0 0 0 0 0 0 0 0 0 0 13
14 17167.2 16772.1 0 0 1 0 0 0 0 0 0 0 0 0 14
15 19630.0 19106.9 0 0 0 1 0 0 0 0 0 0 0 0 15
16 17183.6 16721.3 0 0 0 0 1 0 0 0 0 0 0 0 16
17 18344.7 18161.3 0 0 0 0 0 1 0 0 0 0 0 0 17
18 19301.4 18509.9 0 0 0 0 0 0 1 0 0 0 0 0 18
19 18147.5 17802.7 0 0 0 0 0 0 0 1 0 0 0 0 19
20 16192.9 16409.9 0 0 0 0 0 0 0 0 1 0 0 0 20
21 18374.4 17967.7 0 0 0 0 0 0 0 0 0 1 0 0 21
22 20515.2 20286.6 0 0 0 0 0 0 0 0 0 0 1 0 22
23 18957.2 19537.3 0 0 0 0 0 0 0 0 0 0 0 1 23
24 16471.5 18021.9 0 0 0 0 0 0 0 0 0 0 0 0 24
25 18746.8 20194.3 0 1 0 0 0 0 0 0 0 0 0 0 25
26 19009.5 19049.6 0 0 1 0 0 0 0 0 0 0 0 0 26
27 19211.2 20244.7 0 0 0 1 0 0 0 0 0 0 0 0 27
28 20547.7 21473.3 0 0 0 0 1 0 0 0 0 0 0 0 28
29 19325.8 19673.6 0 0 0 0 0 1 0 0 0 0 0 0 29
30 20605.5 21053.2 0 0 0 0 0 0 1 0 0 0 0 0 30
31 20056.9 20159.5 0 0 0 0 0 0 0 1 0 0 0 0 31
32 16141.4 18203.6 0 0 0 0 0 0 0 0 1 0 0 0 32
33 20359.8 21289.5 0 0 0 0 0 0 0 0 0 1 0 0 33
34 19711.6 20432.3 1 0 0 0 0 0 0 0 0 0 1 0 34
35 15638.6 17180.4 1 0 0 0 0 0 0 0 0 0 0 1 35
36 14384.5 15816.8 1 0 0 0 0 0 0 0 0 0 0 0 36
37 13855.6 15071.8 1 1 0 0 0 0 0 0 0 0 0 0 37
38 14308.3 14521.1 1 0 1 0 0 0 0 0 0 0 0 0 38
39 15290.6 15668.8 1 0 0 1 0 0 0 0 0 0 0 0 39
40 14423.8 14346.9 1 0 0 0 1 0 0 0 0 0 0 0 40
41 13779.7 13881.0 1 0 0 0 0 1 0 0 0 0 0 0 41
42 15686.3 15465.9 1 0 0 0 0 0 1 0 0 0 0 0 42
43 14733.8 14238.2 1 0 0 0 0 0 0 1 0 0 0 0 43
44 12522.5 13557.7 1 0 0 0 0 0 0 0 1 0 0 0 44
45 16189.4 16127.6 1 0 0 0 0 0 0 0 0 1 0 0 45
46 16059.1 16793.9 1 0 0 0 0 0 0 0 0 0 1 0 46
47 16007.1 16014.0 1 0 0 0 0 0 0 0 0 0 0 1 47
48 15806.8 16867.9 1 0 0 0 0 0 0 0 0 0 0 0 48
49 15160.0 16014.6 0 1 0 0 0 0 0 0 0 0 0 0 49
50 15692.1 15878.6 0 0 1 0 0 0 0 0 0 0 0 0 50
51 18908.9 18664.9 0 0 0 1 0 0 0 0 0 0 0 0 51
52 16969.9 17962.5 0 0 0 0 1 0 0 0 0 0 0 0 52
53 16997.5 17332.7 0 0 0 0 0 1 0 0 0 0 0 0 53
54 19858.9 19542.1 0 0 0 0 0 0 1 0 0 0 0 0 54
55 17681.2 17203.6 0 0 0 0 0 0 0 1 0 0 0 0 55
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) invoer crisis M1 M2 M3
3664.9954 0.7581 -755.1942 21.5711 712.1540 1108.9865
M4 M5 M6 M7 M8 M9
658.0625 943.4359 1543.2113 1336.5540 -396.9397 1307.0077
M10 M11 t
1409.1436 881.6974 -9.7033
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-837.95 -232.62 37.52 233.94 678.37
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3664.99536 797.54309 4.595 4.26e-05 ***
invoer 0.75811 0.04339 17.473 < 2e-16 ***
crisis -755.19417 204.57396 -3.692 0.000665 ***
M1 21.57110 284.60683 0.076 0.939962
M2 712.15399 287.24263 2.479 0.017478 *
M3 1108.98646 288.25277 3.847 0.000420 ***
M4 658.06250 284.80601 2.311 0.026096 *
M5 943.43588 285.08977 3.309 0.001987 **
M6 1543.21132 287.19982 5.373 3.58e-06 ***
M7 1336.55398 287.30268 4.652 3.56e-05 ***
M8 -396.93967 304.79381 -1.302 0.200255
M9 1307.00770 300.05155 4.356 8.97e-05 ***
M10 1409.14358 310.61332 4.537 5.11e-05 ***
M11 881.69742 299.41356 2.945 0.005363 **
t -9.70329 4.31070 -2.251 0.029951 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 419.3 on 40 degrees of freedom
Multiple R-squared: 0.9681, Adjusted R-squared: 0.9569
F-statistic: 86.57 on 14 and 40 DF, p-value: < 2.2e-16
> if (n > n25) {
+ kp3 <- k + 3
+ nmkm3 <- n - k - 3
+ gqarr <- array(NA, dim=c(nmkm3-kp3+1,3))
+ numgqtests <- 0
+ numsignificant1 <- 0
+ numsignificant5 <- 0
+ numsignificant10 <- 0
+ for (mypoint in kp3:nmkm3) {
+ j <- 0
+ numgqtests <- numgqtests + 1
+ for (myalt in c('greater', 'two.sided', 'less')) {
+ j <- j + 1
+ gqarr[mypoint-kp3+1,j] <- gqtest(mylm, point=mypoint, alternative=myalt)$p.value
+ }
+ if (gqarr[mypoint-kp3+1,2] < 0.01) numsignificant1 <- numsignificant1 + 1
+ if (gqarr[mypoint-kp3+1,2] < 0.05) numsignificant5 <- numsignificant5 + 1
+ if (gqarr[mypoint-kp3+1,2] < 0.10) numsignificant10 <- numsignificant10 + 1
+ }
+ gqarr
+ }
[,1] [,2] [,3]
[1,] 0.2950825 0.5901651 0.7049175
[2,] 0.1575580 0.3151160 0.8424420
[3,] 0.1550811 0.3101622 0.8449189
[4,] 0.2569608 0.5139216 0.7430392
[5,] 0.2542235 0.5084471 0.7457765
[6,] 0.4472282 0.8944564 0.5527718
[7,] 0.5475664 0.9048672 0.4524336
[8,] 0.4329202 0.8658404 0.5670798
[9,] 0.5064725 0.9870549 0.4935275
[10,] 0.7581450 0.4837100 0.2418550
[11,] 0.7051647 0.5896705 0.2948353
[12,] 0.7066966 0.5866068 0.2933034
[13,] 0.6188481 0.7623037 0.3811519
[14,] 0.5554106 0.8891787 0.4445894
[15,] 0.5407634 0.9184732 0.4592366
[16,] 0.4442002 0.8884003 0.5557998
[17,] 0.6482891 0.7034218 0.3517109
[18,] 0.6409604 0.7180791 0.3590396
[19,] 0.4940153 0.9880306 0.5059847
[20,] 0.3265741 0.6531482 0.6734259
> postscript(file="/var/www/html/freestat/rcomp/tmp/12wy81291119898.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/2vnyt1291119898.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/3vnyt1291119898.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/4vnyt1291119898.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/55efv1291119898.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 = 55
Frequency = 1
1 2 3 4 5 6
-287.171747 -462.493592 -119.106246 -457.696327 218.928128 -238.510898
7 8 9 10 11 12
-188.492897 -220.188143 264.459372 37.524864 484.468181 215.113824
13 14 15 16 17 18
628.715137 210.769843 516.401242 339.180004 132.928921 235.279008
19 20 21 22 23 24
-166.123698 678.371357 -15.359284 275.022679 -177.774696 -623.231358
25 26 27 28 29 30
-6.721234 442.909739 -648.538844 217.172253 83.975966 -272.287212
31 32 33 34 35 36
72.997933 -616.514258 -431.815539 232.599483 -837.947303 -166.885337
37 38 39 40 41 42
-142.859849 46.452716 -228.461368 367.513873 -199.051932 -84.055463
43 44 45 46 47 48
110.539023 158.331044 182.715451 -545.147027 531.253819 575.002870
49 50 51 52 53 54
-191.962308 -237.638706 479.705216 -466.169803 -236.781084 359.574566
55
171.079639
> postscript(file="/var/www/html/freestat/rcomp/tmp/65efv1291119898.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 = 55
Frequency = 1
lag(myerror, k = 1) myerror
0 -287.171747 NA
1 -462.493592 -287.171747
2 -119.106246 -462.493592
3 -457.696327 -119.106246
4 218.928128 -457.696327
5 -238.510898 218.928128
6 -188.492897 -238.510898
7 -220.188143 -188.492897
8 264.459372 -220.188143
9 37.524864 264.459372
10 484.468181 37.524864
11 215.113824 484.468181
12 628.715137 215.113824
13 210.769843 628.715137
14 516.401242 210.769843
15 339.180004 516.401242
16 132.928921 339.180004
17 235.279008 132.928921
18 -166.123698 235.279008
19 678.371357 -166.123698
20 -15.359284 678.371357
21 275.022679 -15.359284
22 -177.774696 275.022679
23 -623.231358 -177.774696
24 -6.721234 -623.231358
25 442.909739 -6.721234
26 -648.538844 442.909739
27 217.172253 -648.538844
28 83.975966 217.172253
29 -272.287212 83.975966
30 72.997933 -272.287212
31 -616.514258 72.997933
32 -431.815539 -616.514258
33 232.599483 -431.815539
34 -837.947303 232.599483
35 -166.885337 -837.947303
36 -142.859849 -166.885337
37 46.452716 -142.859849
38 -228.461368 46.452716
39 367.513873 -228.461368
40 -199.051932 367.513873
41 -84.055463 -199.051932
42 110.539023 -84.055463
43 158.331044 110.539023
44 182.715451 158.331044
45 -545.147027 182.715451
46 531.253819 -545.147027
47 575.002870 531.253819
48 -191.962308 575.002870
49 -237.638706 -191.962308
50 479.705216 -237.638706
51 -466.169803 479.705216
52 -236.781084 -466.169803
53 359.574566 -236.781084
54 171.079639 359.574566
55 NA 171.079639
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -462.493592 -287.171747
[2,] -119.106246 -462.493592
[3,] -457.696327 -119.106246
[4,] 218.928128 -457.696327
[5,] -238.510898 218.928128
[6,] -188.492897 -238.510898
[7,] -220.188143 -188.492897
[8,] 264.459372 -220.188143
[9,] 37.524864 264.459372
[10,] 484.468181 37.524864
[11,] 215.113824 484.468181
[12,] 628.715137 215.113824
[13,] 210.769843 628.715137
[14,] 516.401242 210.769843
[15,] 339.180004 516.401242
[16,] 132.928921 339.180004
[17,] 235.279008 132.928921
[18,] -166.123698 235.279008
[19,] 678.371357 -166.123698
[20,] -15.359284 678.371357
[21,] 275.022679 -15.359284
[22,] -177.774696 275.022679
[23,] -623.231358 -177.774696
[24,] -6.721234 -623.231358
[25,] 442.909739 -6.721234
[26,] -648.538844 442.909739
[27,] 217.172253 -648.538844
[28,] 83.975966 217.172253
[29,] -272.287212 83.975966
[30,] 72.997933 -272.287212
[31,] -616.514258 72.997933
[32,] -431.815539 -616.514258
[33,] 232.599483 -431.815539
[34,] -837.947303 232.599483
[35,] -166.885337 -837.947303
[36,] -142.859849 -166.885337
[37,] 46.452716 -142.859849
[38,] -228.461368 46.452716
[39,] 367.513873 -228.461368
[40,] -199.051932 367.513873
[41,] -84.055463 -199.051932
[42,] 110.539023 -84.055463
[43,] 158.331044 110.539023
[44,] 182.715451 158.331044
[45,] -545.147027 182.715451
[46,] 531.253819 -545.147027
[47,] 575.002870 531.253819
[48,] -191.962308 575.002870
[49,] -237.638706 -191.962308
[50,] 479.705216 -237.638706
[51,] -466.169803 479.705216
[52,] -236.781084 -466.169803
[53,] 359.574566 -236.781084
[54,] 171.079639 359.574566
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -462.493592 -287.171747
2 -119.106246 -462.493592
3 -457.696327 -119.106246
4 218.928128 -457.696327
5 -238.510898 218.928128
6 -188.492897 -238.510898
7 -220.188143 -188.492897
8 264.459372 -220.188143
9 37.524864 264.459372
10 484.468181 37.524864
11 215.113824 484.468181
12 628.715137 215.113824
13 210.769843 628.715137
14 516.401242 210.769843
15 339.180004 516.401242
16 132.928921 339.180004
17 235.279008 132.928921
18 -166.123698 235.279008
19 678.371357 -166.123698
20 -15.359284 678.371357
21 275.022679 -15.359284
22 -177.774696 275.022679
23 -623.231358 -177.774696
24 -6.721234 -623.231358
25 442.909739 -6.721234
26 -648.538844 442.909739
27 217.172253 -648.538844
28 83.975966 217.172253
29 -272.287212 83.975966
30 72.997933 -272.287212
31 -616.514258 72.997933
32 -431.815539 -616.514258
33 232.599483 -431.815539
34 -837.947303 232.599483
35 -166.885337 -837.947303
36 -142.859849 -166.885337
37 46.452716 -142.859849
38 -228.461368 46.452716
39 367.513873 -228.461368
40 -199.051932 367.513873
41 -84.055463 -199.051932
42 110.539023 -84.055463
43 158.331044 110.539023
44 182.715451 158.331044
45 -545.147027 182.715451
46 531.253819 -545.147027
47 575.002870 531.253819
48 -191.962308 575.002870
49 -237.638706 -191.962308
50 479.705216 -237.638706
51 -466.169803 479.705216
52 -236.781084 -466.169803
53 359.574566 -236.781084
54 171.079639 359.574566
> 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/7yoey1291119898.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/8yoey1291119898.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/99xdj1291119898.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/109xdj1291119898.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/11ufu71291119898.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/12ggad1291119898.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/134hpp1291119898.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/14x8pa1291119898.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/151r5x1291119898.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/16xi3o1291119898.tab")
+ }
>
> try(system("convert tmp/12wy81291119898.ps tmp/12wy81291119898.png",intern=TRUE))
character(0)
> try(system("convert tmp/2vnyt1291119898.ps tmp/2vnyt1291119898.png",intern=TRUE))
character(0)
> try(system("convert tmp/3vnyt1291119898.ps tmp/3vnyt1291119898.png",intern=TRUE))
character(0)
> try(system("convert tmp/4vnyt1291119898.ps tmp/4vnyt1291119898.png",intern=TRUE))
character(0)
> try(system("convert tmp/55efv1291119898.ps tmp/55efv1291119898.png",intern=TRUE))
character(0)
> try(system("convert tmp/65efv1291119898.ps tmp/65efv1291119898.png",intern=TRUE))
character(0)
> try(system("convert tmp/7yoey1291119898.ps tmp/7yoey1291119898.png",intern=TRUE))
character(0)
> try(system("convert tmp/8yoey1291119898.ps tmp/8yoey1291119898.png",intern=TRUE))
character(0)
> try(system("convert tmp/99xdj1291119898.ps tmp/99xdj1291119898.png",intern=TRUE))
character(0)
> try(system("convert tmp/109xdj1291119898.ps tmp/109xdj1291119898.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.668 2.441 4.014