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(235.1,9700,280.7,9081,264.6,9084,240.7,9743,201.4,8587,240.8,9731,241.1,9563,223.8,9998,206.1,9437,174.7,10038,203.3,9918,220.5,9252,299.5,9737,347.4,9035,338.3,9133,327.7,9487,351.6,8700,396.6,9627,438.8,8947,395.6,9283,363.5,8829,378.8,9947,357,9628,369,9318,464.8,9605,479.1,8640,431.3,9214,366.5,9567,326.3,8547,355.1,9185,331.6,9470,261.3,9123,249,9278,205.5,10170,235.6,9434,240.9,9655,264.9,9429,253.8,8739,232.3,9552,193.8,9687,177,9019,213.2,9672,207.2,9206,180.6,9069,188.6,9788,175.4,10312,199,10105,179.6,9863,225.8,9656,234,9295,200.2,9946,183.6,9701,178.2,9049,203.2,10190,208.5,9706,191.8,9765,172.8,9893,148,9994,159.4,10433,154.5,10073,213.2,10112,196.4,9266,182.8,9820,176.4,10097,153.6,9115,173.2,10411,171,9678,151.2,10408,161.9,10153,157.2,10368,201.7,10581,236.4,10597,356.1,10680,398.3,9738,403.7,9556),dim=c(2,75),dimnames=list(c('unemployment','birth'),1:75))
> y <- array(NA,dim=c(2,75),dimnames=list(c('unemployment','birth'),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 = 'Do not include Seasonal Dummies'
> par1 = '1'
> #'GNU S' R Code compiled by R2WASP v. 1.0.44 ()
> #Author: Prof. Dr. P. Wessa
> #To cite this work: AUTHOR(S), (YEAR), YOUR SOFTWARE TITLE (vNUMBER) in Free Statistics Software (v$_version), Office for Research Development and Education, URL http://www.wessa.net/rwasp_YOURPAGE.wasp/
> #Source of accompanying publication: Office for Research, Development, and Education
> #Technical description: Write here your technical program description (don't use hard returns!)
> library(lattice)
> library(lmtest)
Loading required package: zoo
Attaching package: 'zoo'
The following object(s) are masked from package:base :
as.Date.numeric
> n25 <- 25 #minimum number of obs. for Goldfeld-Quandt test
> par1 <- as.numeric(par1)
> x <- t(y)
> k <- length(x[1,])
> n <- length(x[,1])
> x1 <- cbind(x[,par1], x[,1:k!=par1])
> mycolnames <- c(colnames(x)[par1], colnames(x)[1:k!=par1])
> colnames(x1) <- mycolnames #colnames(x)[par1]
> x <- x1
> if (par3 == 'First Differences'){
+ x2 <- array(0, dim=c(n-1,k), dimnames=list(1:(n-1), paste('(1-B)',colnames(x),sep='')))
+ for (i in 1:n-1) {
+ for (j in 1:k) {
+ x2[i,j] <- x[i+1,j] - x[i,j]
+ }
+ }
+ x <- x2
+ }
> if (par2 == 'Include Monthly Dummies'){
+ x2 <- array(0, dim=c(n,11), dimnames=list(1:n, paste('M', seq(1:11), sep ='')))
+ for (i in 1:11){
+ x2[seq(i,n,12),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> if (par2 == 'Include Quarterly Dummies'){
+ x2 <- array(0, dim=c(n,3), dimnames=list(1:n, paste('Q', seq(1:3), sep ='')))
+ for (i in 1:3){
+ x2[seq(i,n,4),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> k <- length(x[1,])
> if (par3 == 'Linear Trend'){
+ x <- cbind(x, c(1:n))
+ colnames(x)[k+1] <- 't'
+ }
> x
unemployment birth
1 235.1 9700
2 280.7 9081
3 264.6 9084
4 240.7 9743
5 201.4 8587
6 240.8 9731
7 241.1 9563
8 223.8 9998
9 206.1 9437
10 174.7 10038
11 203.3 9918
12 220.5 9252
13 299.5 9737
14 347.4 9035
15 338.3 9133
16 327.7 9487
17 351.6 8700
18 396.6 9627
19 438.8 8947
20 395.6 9283
21 363.5 8829
22 378.8 9947
23 357.0 9628
24 369.0 9318
25 464.8 9605
26 479.1 8640
27 431.3 9214
28 366.5 9567
29 326.3 8547
30 355.1 9185
31 331.6 9470
32 261.3 9123
33 249.0 9278
34 205.5 10170
35 235.6 9434
36 240.9 9655
37 264.9 9429
38 253.8 8739
39 232.3 9552
40 193.8 9687
41 177.0 9019
42 213.2 9672
43 207.2 9206
44 180.6 9069
45 188.6 9788
46 175.4 10312
47 199.0 10105
48 179.6 9863
49 225.8 9656
50 234.0 9295
51 200.2 9946
52 183.6 9701
53 178.2 9049
54 203.2 10190
55 208.5 9706
56 191.8 9765
57 172.8 9893
58 148.0 9994
59 159.4 10433
60 154.5 10073
61 213.2 10112
62 196.4 9266
63 182.8 9820
64 176.4 10097
65 153.6 9115
66 173.2 10411
67 171.0 9678
68 151.2 10408
69 161.9 10153
70 157.2 10368
71 201.7 10581
72 236.4 10597
73 356.1 10680
74 398.3 9738
75 403.7 9556
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) birth
943.4030 -0.0717
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-136.26 -53.27 -25.27 52.30 210.07
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 943.40303 177.73912 5.308 1.15e-06 ***
birth -0.07170 0.01848 -3.880 0.000227 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 79.97 on 73 degrees of freedom
Multiple R-squared: 0.171, Adjusted R-squared: 0.1596
F-statistic: 15.05 on 1 and 73 DF, p-value: 0.0002267
> 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.0882557727 0.1765115455 0.91174423
[2,] 0.0304899669 0.0609799339 0.96951003
[3,] 0.0093999963 0.0187999926 0.99060000
[4,] 0.0036120603 0.0072241206 0.99638794
[5,] 0.0021963512 0.0043927024 0.99780365
[6,] 0.0032421952 0.0064843904 0.99675780
[7,] 0.0013031469 0.0026062939 0.99869685
[8,] 0.0004902089 0.0009804178 0.99950979
[9,] 0.0015745562 0.0031491124 0.99842544
[10,] 0.0080084817 0.0160169634 0.99199152
[11,] 0.0120976776 0.0241953553 0.98790232
[12,] 0.0157584251 0.0315168502 0.98424157
[13,] 0.0130870861 0.0261741722 0.98691291
[14,] 0.0697735838 0.1395471676 0.93022642
[15,] 0.1684625334 0.3369250668 0.83153747
[16,] 0.2283503217 0.4567006433 0.77164968
[17,] 0.1889908261 0.3779816522 0.81100917
[18,] 0.3148611981 0.6297223963 0.68513880
[19,] 0.3342005850 0.6684011700 0.66579942
[20,] 0.3404199369 0.6808398738 0.65958006
[21,] 0.6967834226 0.6064331547 0.30321658
[22,] 0.8553861612 0.2892276776 0.14461384
[23,] 0.9406597036 0.1186805927 0.05934030
[24,] 0.9618260768 0.0763478463 0.03817392
[25,] 0.9579593946 0.0840812109 0.04204061
[26,] 0.9690239117 0.0619521766 0.03097609
[27,] 0.9755476032 0.0489047937 0.02445240
[28,] 0.9722379025 0.0555241951 0.02776210
[29,] 0.9674850535 0.0650298930 0.03251495
[30,] 0.9569505045 0.0860989911 0.04304950
[31,] 0.9483849186 0.1032301627 0.05161508
[32,] 0.9348817725 0.1302364549 0.06511823
[33,] 0.9239594247 0.1520811506 0.07604058
[34,] 0.9238201054 0.1523597891 0.07617989
[35,] 0.9075131319 0.1849737362 0.09248687
[36,] 0.8946505501 0.2106988999 0.10534945
[37,] 0.9111483465 0.1777033070 0.08885165
[38,] 0.8892356305 0.2215287390 0.11076437
[39,] 0.8755391976 0.2489216047 0.12446080
[40,] 0.8782587772 0.2434824456 0.12174122
[41,] 0.8539317863 0.2921364274 0.14606821
[42,] 0.8210414023 0.3579171953 0.17895860
[43,] 0.7748787710 0.4502424579 0.22512123
[44,] 0.7385247699 0.5229504601 0.26147523
[45,] 0.6836416025 0.6327167950 0.31635840
[46,] 0.6329051993 0.7341896015 0.36709480
[47,] 0.5674331769 0.8651336461 0.43256682
[48,] 0.5166947942 0.9666104116 0.48330521
[49,] 0.5054349233 0.9891301533 0.49456508
[50,] 0.4304134336 0.8608268671 0.56958657
[51,] 0.3622514155 0.7245028310 0.63774858
[52,] 0.3029591322 0.6059182645 0.69704087
[53,] 0.2571796961 0.5143593921 0.74282030
[54,] 0.2364276875 0.4728553751 0.76357231
[55,] 0.1951917149 0.3903834298 0.80480829
[56,] 0.1709413110 0.3418826219 0.82905869
[57,] 0.1221072222 0.2442144444 0.87789278
[58,] 0.0961039537 0.1922079073 0.90389605
[59,] 0.0720174616 0.1440349231 0.92798254
[60,] 0.0523531231 0.1047062462 0.94764688
[61,] 0.0999267752 0.1998535504 0.90007322
[62,] 0.0707903087 0.1415806174 0.92920969
[63,] 0.1327459114 0.2654918229 0.86725409
[64,] 0.1290769899 0.2581539799 0.87092301
[65,] 0.2248577840 0.4497155679 0.77514222
[66,] 0.4236953212 0.8473906424 0.57630468
> postscript(file="/var/www/html/freestat/rcomp/tmp/13fsh1293205404.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/23fsh1293205404.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/3eoa21293205404.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/4eoa21293205404.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/5eoa21293205404.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
-12.819558 -11.601442 -27.486344 -4.136487 -126.320909 -4.896879
7 8 9 10 11 12
-16.642366 -2.753159 -60.676481 -48.985186 -28.989105 -59.540857
13 14 15 16 17 18
54.233317 51.800389 49.726923 64.508485 31.981115 143.446391
19 20 21 22 23 24
136.890848 117.781822 53.130328 148.590176 103.918090 93.691299
25 26 27 28 29 30
210.069006 155.179155 148.534569 109.044431 -4.288882 70.255288
31 32 33 34 35 36
67.189597 -27.990070 -29.176674 -8.720875 -31.391579 -10.246028
37 38 39 40 41 42
-2.450076 -63.022612 -26.231059 -55.051650 -119.746800 -36.727139
43 44 45 46 47 48
-76.139026 -112.561834 -53.010017 -28.639570 -19.881331 -56.632568
49 50 51 52 53 54
-25.274329 -42.957786 -30.081524 -64.247859 -116.395820 -9.586888
55 56 57 58 59 60
-38.989362 -51.459102 -61.281588 -78.839956 -35.963952 -66.675709
61 62 63 64 65 66
-5.179436 -82.637066 -56.515639 -43.054925 -136.263665 -23.741337
67 68 69 70 71 72
-78.496943 -45.956435 -53.539763 -42.824408 16.947549 52.794738
73 74 75
178.445782 153.105016 145.455739
> postscript(file="/var/www/html/freestat/rcomp/tmp/66yr51293205404.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 -12.819558 NA
1 -11.601442 -12.819558
2 -27.486344 -11.601442
3 -4.136487 -27.486344
4 -126.320909 -4.136487
5 -4.896879 -126.320909
6 -16.642366 -4.896879
7 -2.753159 -16.642366
8 -60.676481 -2.753159
9 -48.985186 -60.676481
10 -28.989105 -48.985186
11 -59.540857 -28.989105
12 54.233317 -59.540857
13 51.800389 54.233317
14 49.726923 51.800389
15 64.508485 49.726923
16 31.981115 64.508485
17 143.446391 31.981115
18 136.890848 143.446391
19 117.781822 136.890848
20 53.130328 117.781822
21 148.590176 53.130328
22 103.918090 148.590176
23 93.691299 103.918090
24 210.069006 93.691299
25 155.179155 210.069006
26 148.534569 155.179155
27 109.044431 148.534569
28 -4.288882 109.044431
29 70.255288 -4.288882
30 67.189597 70.255288
31 -27.990070 67.189597
32 -29.176674 -27.990070
33 -8.720875 -29.176674
34 -31.391579 -8.720875
35 -10.246028 -31.391579
36 -2.450076 -10.246028
37 -63.022612 -2.450076
38 -26.231059 -63.022612
39 -55.051650 -26.231059
40 -119.746800 -55.051650
41 -36.727139 -119.746800
42 -76.139026 -36.727139
43 -112.561834 -76.139026
44 -53.010017 -112.561834
45 -28.639570 -53.010017
46 -19.881331 -28.639570
47 -56.632568 -19.881331
48 -25.274329 -56.632568
49 -42.957786 -25.274329
50 -30.081524 -42.957786
51 -64.247859 -30.081524
52 -116.395820 -64.247859
53 -9.586888 -116.395820
54 -38.989362 -9.586888
55 -51.459102 -38.989362
56 -61.281588 -51.459102
57 -78.839956 -61.281588
58 -35.963952 -78.839956
59 -66.675709 -35.963952
60 -5.179436 -66.675709
61 -82.637066 -5.179436
62 -56.515639 -82.637066
63 -43.054925 -56.515639
64 -136.263665 -43.054925
65 -23.741337 -136.263665
66 -78.496943 -23.741337
67 -45.956435 -78.496943
68 -53.539763 -45.956435
69 -42.824408 -53.539763
70 16.947549 -42.824408
71 52.794738 16.947549
72 178.445782 52.794738
73 153.105016 178.445782
74 145.455739 153.105016
75 NA 145.455739
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -11.601442 -12.819558
[2,] -27.486344 -11.601442
[3,] -4.136487 -27.486344
[4,] -126.320909 -4.136487
[5,] -4.896879 -126.320909
[6,] -16.642366 -4.896879
[7,] -2.753159 -16.642366
[8,] -60.676481 -2.753159
[9,] -48.985186 -60.676481
[10,] -28.989105 -48.985186
[11,] -59.540857 -28.989105
[12,] 54.233317 -59.540857
[13,] 51.800389 54.233317
[14,] 49.726923 51.800389
[15,] 64.508485 49.726923
[16,] 31.981115 64.508485
[17,] 143.446391 31.981115
[18,] 136.890848 143.446391
[19,] 117.781822 136.890848
[20,] 53.130328 117.781822
[21,] 148.590176 53.130328
[22,] 103.918090 148.590176
[23,] 93.691299 103.918090
[24,] 210.069006 93.691299
[25,] 155.179155 210.069006
[26,] 148.534569 155.179155
[27,] 109.044431 148.534569
[28,] -4.288882 109.044431
[29,] 70.255288 -4.288882
[30,] 67.189597 70.255288
[31,] -27.990070 67.189597
[32,] -29.176674 -27.990070
[33,] -8.720875 -29.176674
[34,] -31.391579 -8.720875
[35,] -10.246028 -31.391579
[36,] -2.450076 -10.246028
[37,] -63.022612 -2.450076
[38,] -26.231059 -63.022612
[39,] -55.051650 -26.231059
[40,] -119.746800 -55.051650
[41,] -36.727139 -119.746800
[42,] -76.139026 -36.727139
[43,] -112.561834 -76.139026
[44,] -53.010017 -112.561834
[45,] -28.639570 -53.010017
[46,] -19.881331 -28.639570
[47,] -56.632568 -19.881331
[48,] -25.274329 -56.632568
[49,] -42.957786 -25.274329
[50,] -30.081524 -42.957786
[51,] -64.247859 -30.081524
[52,] -116.395820 -64.247859
[53,] -9.586888 -116.395820
[54,] -38.989362 -9.586888
[55,] -51.459102 -38.989362
[56,] -61.281588 -51.459102
[57,] -78.839956 -61.281588
[58,] -35.963952 -78.839956
[59,] -66.675709 -35.963952
[60,] -5.179436 -66.675709
[61,] -82.637066 -5.179436
[62,] -56.515639 -82.637066
[63,] -43.054925 -56.515639
[64,] -136.263665 -43.054925
[65,] -23.741337 -136.263665
[66,] -78.496943 -23.741337
[67,] -45.956435 -78.496943
[68,] -53.539763 -45.956435
[69,] -42.824408 -53.539763
[70,] 16.947549 -42.824408
[71,] 52.794738 16.947549
[72,] 178.445782 52.794738
[73,] 153.105016 178.445782
[74,] 145.455739 153.105016
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -11.601442 -12.819558
2 -27.486344 -11.601442
3 -4.136487 -27.486344
4 -126.320909 -4.136487
5 -4.896879 -126.320909
6 -16.642366 -4.896879
7 -2.753159 -16.642366
8 -60.676481 -2.753159
9 -48.985186 -60.676481
10 -28.989105 -48.985186
11 -59.540857 -28.989105
12 54.233317 -59.540857
13 51.800389 54.233317
14 49.726923 51.800389
15 64.508485 49.726923
16 31.981115 64.508485
17 143.446391 31.981115
18 136.890848 143.446391
19 117.781822 136.890848
20 53.130328 117.781822
21 148.590176 53.130328
22 103.918090 148.590176
23 93.691299 103.918090
24 210.069006 93.691299
25 155.179155 210.069006
26 148.534569 155.179155
27 109.044431 148.534569
28 -4.288882 109.044431
29 70.255288 -4.288882
30 67.189597 70.255288
31 -27.990070 67.189597
32 -29.176674 -27.990070
33 -8.720875 -29.176674
34 -31.391579 -8.720875
35 -10.246028 -31.391579
36 -2.450076 -10.246028
37 -63.022612 -2.450076
38 -26.231059 -63.022612
39 -55.051650 -26.231059
40 -119.746800 -55.051650
41 -36.727139 -119.746800
42 -76.139026 -36.727139
43 -112.561834 -76.139026
44 -53.010017 -112.561834
45 -28.639570 -53.010017
46 -19.881331 -28.639570
47 -56.632568 -19.881331
48 -25.274329 -56.632568
49 -42.957786 -25.274329
50 -30.081524 -42.957786
51 -64.247859 -30.081524
52 -116.395820 -64.247859
53 -9.586888 -116.395820
54 -38.989362 -9.586888
55 -51.459102 -38.989362
56 -61.281588 -51.459102
57 -78.839956 -61.281588
58 -35.963952 -78.839956
59 -66.675709 -35.963952
60 -5.179436 -66.675709
61 -82.637066 -5.179436
62 -56.515639 -82.637066
63 -43.054925 -56.515639
64 -136.263665 -43.054925
65 -23.741337 -136.263665
66 -78.496943 -23.741337
67 -45.956435 -78.496943
68 -53.539763 -45.956435
69 -42.824408 -53.539763
70 16.947549 -42.824408
71 52.794738 16.947549
72 178.445782 52.794738
73 153.105016 178.445782
74 145.455739 153.105016
> 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/74afe1293205404.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/84afe1293205404.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/9ay8t1293205404.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/10ay8t1293205404.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/11dzoh1293205404.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/12yh441293205404.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/13niky1293205404.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/14yrj11293205404.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/151azp1293205404.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/16nbgv1293205404.tab")
+ }
>
> try(system("convert tmp/13fsh1293205404.ps tmp/13fsh1293205404.png",intern=TRUE))
character(0)
> try(system("convert tmp/23fsh1293205404.ps tmp/23fsh1293205404.png",intern=TRUE))
character(0)
> try(system("convert tmp/3eoa21293205404.ps tmp/3eoa21293205404.png",intern=TRUE))
character(0)
> try(system("convert tmp/4eoa21293205404.ps tmp/4eoa21293205404.png",intern=TRUE))
character(0)
> try(system("convert tmp/5eoa21293205404.ps tmp/5eoa21293205404.png",intern=TRUE))
character(0)
> try(system("convert tmp/66yr51293205404.ps tmp/66yr51293205404.png",intern=TRUE))
character(0)
> try(system("convert tmp/74afe1293205404.ps tmp/74afe1293205404.png",intern=TRUE))
character(0)
> try(system("convert tmp/84afe1293205404.ps tmp/84afe1293205404.png",intern=TRUE))
character(0)
> try(system("convert tmp/9ay8t1293205404.ps tmp/9ay8t1293205404.png",intern=TRUE))
character(0)
> try(system("convert tmp/10ay8t1293205404.ps tmp/10ay8t1293205404.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
4.082 2.554 4.466