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(14731798.37,0,16471559.62,0,15213975.95,0,17637387.4,0,17972385.83,0,16896235.55,0,16697955.94,0,19691579.52,0,15930700.75,0,17444615.98,0,17699369.88,0,15189796.81,0,15672722.75,0,17180794.3,0,17664893.45,0,17862884.98,0,16162288.88,0,17463628.82,0,16772112.17,0,19106861.48,0,16721314.25,0,18161267.85,0,18509941.2,0,17802737.97,0,16409869.75,0,17967742.04,0,20286602.27,0,19537280.81,0,18021889.62,0,20194317.23,0,19049596.62,0,20244720.94,0,21473302.24,0,19673603.19,0,21053177.29,0,20159479.84,0,18203628.31,0,21289464.94,0,20432335.71,1,17180395.07,1,15816786.32,1,15071819.75,1,14521120.61,1,15668789.39,1,14346884.11,1,13881008.13,1,15465943.69,1,14238232.92,1,13557713.21,1,16127590.29,1,16793894.2,1,16014007.43,1,16867867.15,1,16014583.21,1,15878594.85,1,18664899.14,1,17962530.06,1,17332692.2,1,19542066.35,1,17203555.19,1),dim=c(2,60),dimnames=list(c('Y','X'),1:60))
> y <- array(NA,dim=c(2,60),dimnames=list(c('Y','X'),1:60))
> 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
Y X M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 14731798 0 1 0 0 0 0 0 0 0 0 0 0
2 16471560 0 0 1 0 0 0 0 0 0 0 0 0
3 15213976 0 0 0 1 0 0 0 0 0 0 0 0
4 17637387 0 0 0 0 1 0 0 0 0 0 0 0
5 17972386 0 0 0 0 0 1 0 0 0 0 0 0
6 16896236 0 0 0 0 0 0 1 0 0 0 0 0
7 16697956 0 0 0 0 0 0 0 1 0 0 0 0
8 19691580 0 0 0 0 0 0 0 0 1 0 0 0
9 15930701 0 0 0 0 0 0 0 0 0 1 0 0
10 17444616 0 0 0 0 0 0 0 0 0 0 1 0
11 17699370 0 0 0 0 0 0 0 0 0 0 0 1
12 15189797 0 0 0 0 0 0 0 0 0 0 0 0
13 15672723 0 1 0 0 0 0 0 0 0 0 0 0
14 17180794 0 0 1 0 0 0 0 0 0 0 0 0
15 17664893 0 0 0 1 0 0 0 0 0 0 0 0
16 17862885 0 0 0 0 1 0 0 0 0 0 0 0
17 16162289 0 0 0 0 0 1 0 0 0 0 0 0
18 17463629 0 0 0 0 0 0 1 0 0 0 0 0
19 16772112 0 0 0 0 0 0 0 1 0 0 0 0
20 19106861 0 0 0 0 0 0 0 0 1 0 0 0
21 16721314 0 0 0 0 0 0 0 0 0 1 0 0
22 18161268 0 0 0 0 0 0 0 0 0 0 1 0
23 18509941 0 0 0 0 0 0 0 0 0 0 0 1
24 17802738 0 0 0 0 0 0 0 0 0 0 0 0
25 16409870 0 1 0 0 0 0 0 0 0 0 0 0
26 17967742 0 0 1 0 0 0 0 0 0 0 0 0
27 20286602 0 0 0 1 0 0 0 0 0 0 0 0
28 19537281 0 0 0 0 1 0 0 0 0 0 0 0
29 18021890 0 0 0 0 0 1 0 0 0 0 0 0
30 20194317 0 0 0 0 0 0 1 0 0 0 0 0
31 19049597 0 0 0 0 0 0 0 1 0 0 0 0
32 20244721 0 0 0 0 0 0 0 0 1 0 0 0
33 21473302 0 0 0 0 0 0 0 0 0 1 0 0
34 19673603 0 0 0 0 0 0 0 0 0 0 1 0
35 21053177 0 0 0 0 0 0 0 0 0 0 0 1
36 20159480 0 0 0 0 0 0 0 0 0 0 0 0
37 18203628 0 1 0 0 0 0 0 0 0 0 0 0
38 21289465 0 0 1 0 0 0 0 0 0 0 0 0
39 20432336 1 0 0 1 0 0 0 0 0 0 0 0
40 17180395 1 0 0 0 1 0 0 0 0 0 0 0
41 15816786 1 0 0 0 0 1 0 0 0 0 0 0
42 15071820 1 0 0 0 0 0 1 0 0 0 0 0
43 14521121 1 0 0 0 0 0 0 1 0 0 0 0
44 15668789 1 0 0 0 0 0 0 0 1 0 0 0
45 14346884 1 0 0 0 0 0 0 0 0 1 0 0
46 13881008 1 0 0 0 0 0 0 0 0 0 1 0
47 15465944 1 0 0 0 0 0 0 0 0 0 0 1
48 14238233 1 0 0 0 0 0 0 0 0 0 0 0
49 13557713 1 1 0 0 0 0 0 0 0 0 0 0
50 16127590 1 0 1 0 0 0 0 0 0 0 0 0
51 16793894 1 0 0 1 0 0 0 0 0 0 0 0
52 16014007 1 0 0 0 1 0 0 0 0 0 0 0
53 16867867 1 0 0 0 0 1 0 0 0 0 0 0
54 16014583 1 0 0 0 0 0 1 0 0 0 0 0
55 15878595 1 0 0 0 0 0 0 1 0 0 0 0
56 18664899 1 0 0 0 0 0 0 0 1 0 0 0
57 17962530 1 0 0 0 0 0 0 0 0 1 0 0
58 17332692 1 0 0 0 0 0 0 0 0 0 1 0
59 19542066 1 0 0 0 0 0 0 0 0 0 0 1
60 17203555 1 0 0 0 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
17654550 -1839473 -1571509 520775 1159580 727631
M5 M6 M7 M8 M9 M10
49483 209356 -334885 1756610 368186 379877
M11
1535339
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-3600154 -974386 -207922 1104344 3457679
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 17654550 791301 22.311 < 2e-16 ***
X -1839473 466278 -3.945 0.000265 ***
M1 -1571509 1091530 -1.440 0.156572
M2 520775 1091530 0.477 0.635498
M3 1159580 1087539 1.066 0.291761
M4 727631 1087539 0.669 0.506729
M5 49483 1087539 0.045 0.963902
M6 209356 1087539 0.193 0.848176
M7 -334884 1087539 -0.308 0.759497
M8 1756610 1087539 1.615 0.112958
M9 368186 1087539 0.339 0.736457
M10 379877 1087539 0.349 0.728425
M11 1535339 1087539 1.412 0.164608
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 1720000 on 47 degrees of freedom
Multiple R-squared: 0.3745, Adjusted R-squared: 0.2148
F-statistic: 2.345 on 12 and 47 DF, p-value: 0.01877
> 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.30133106 0.60266213 0.6986689
[2,] 0.28902424 0.57804849 0.7109758
[3,] 0.18077102 0.36154203 0.8192290
[4,] 0.10332543 0.20665085 0.8966746
[5,] 0.05804497 0.11608993 0.9419550
[6,] 0.04651925 0.09303849 0.9534808
[7,] 0.02678035 0.05356069 0.9732197
[8,] 0.01964071 0.03928143 0.9803593
[9,] 0.04991970 0.09983940 0.9500803
[10,] 0.04084394 0.08168789 0.9591561
[11,] 0.04123228 0.08246456 0.9587677
[12,] 0.22320185 0.44640371 0.7767981
[13,] 0.21035579 0.42071157 0.7896442
[14,] 0.20447525 0.40895050 0.7955248
[15,] 0.26558660 0.53117321 0.7344134
[16,] 0.25920677 0.51841354 0.7407932
[17,] 0.20889191 0.41778382 0.7911081
[18,] 0.42168234 0.84336467 0.5783177
[19,] 0.35909353 0.71818705 0.6409065
[20,] 0.35027770 0.70055541 0.6497223
[21,] 0.35496031 0.70992062 0.6450397
[22,] 0.30826025 0.61652050 0.6917398
[23,] 0.31409894 0.62819789 0.6859011
[24,] 0.33289383 0.66578767 0.6671062
[25,] 0.28001742 0.56003483 0.7199826
[26,] 0.20759573 0.41519145 0.7924043
[27,] 0.15183555 0.30367110 0.8481644
[28,] 0.09904816 0.19809632 0.9009518
[29,] 0.08659893 0.17319786 0.9134011
> postscript(file="/var/www/html/freestat/rcomp/tmp/1jnuh1291458480.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/2ufb21291458480.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/3ufb21291458480.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/4ufb21291458480.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/5mosn1291458480.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 = 60
Frequency = 1
1 2 3 4 5 6
-1351242.693 -1703765.203 -3600153.537 -744792.909 268353.099 -967670.533
7 8 9 10 11 12
-621709.269 280420.255 -2092034.703 -589810.661 -1490518.973 -2464752.907
13 14 15 16 17 18
-410318.313 -994530.523 -1149236.037 -519295.329 -1541743.851 -400277.263
19 20 21 22 23 24
-547553.039 -304297.785 -1301421.203 126841.209 -679947.653 148188.253
25 26 27 28 29 30
326828.687 -207582.783 1472472.783 1155100.501 317856.889 2330411.147
31 32 33 34 35 36
1729931.411 833561.675 3450566.787 1639176.549 1863288.437 2504930.123
37 38 39 40 41 42
2120587.247 3114140.117 3457679.150 637687.688 -47773.484 -952613.406
43 44 45 46 47 48
-959071.672 -1902896.948 -1836378.416 -2313945.584 -1884472.236 -1576843.870
49 50 51 52 53 54
-685854.926 -208261.606 -180762.360 -528699.952 1003307.346 -9849.946
55 56 57 58 59 60
398402.568 1093212.802 1779267.534 1137738.486 2191650.424 1388478.400
> postscript(file="/var/www/html/freestat/rcomp/tmp/6mosn1291458480.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 = 60
Frequency = 1
lag(myerror, k = 1) myerror
0 -1351242.693 NA
1 -1703765.203 -1351242.693
2 -3600153.537 -1703765.203
3 -744792.909 -3600153.537
4 268353.099 -744792.909
5 -967670.533 268353.099
6 -621709.269 -967670.533
7 280420.255 -621709.269
8 -2092034.703 280420.255
9 -589810.661 -2092034.703
10 -1490518.973 -589810.661
11 -2464752.907 -1490518.973
12 -410318.313 -2464752.907
13 -994530.523 -410318.313
14 -1149236.037 -994530.523
15 -519295.329 -1149236.037
16 -1541743.851 -519295.329
17 -400277.263 -1541743.851
18 -547553.039 -400277.263
19 -304297.785 -547553.039
20 -1301421.203 -304297.785
21 126841.209 -1301421.203
22 -679947.653 126841.209
23 148188.253 -679947.653
24 326828.687 148188.253
25 -207582.783 326828.687
26 1472472.783 -207582.783
27 1155100.501 1472472.783
28 317856.889 1155100.501
29 2330411.147 317856.889
30 1729931.411 2330411.147
31 833561.675 1729931.411
32 3450566.787 833561.675
33 1639176.549 3450566.787
34 1863288.437 1639176.549
35 2504930.123 1863288.437
36 2120587.247 2504930.123
37 3114140.117 2120587.247
38 3457679.150 3114140.117
39 637687.688 3457679.150
40 -47773.484 637687.688
41 -952613.406 -47773.484
42 -959071.672 -952613.406
43 -1902896.948 -959071.672
44 -1836378.416 -1902896.948
45 -2313945.584 -1836378.416
46 -1884472.236 -2313945.584
47 -1576843.870 -1884472.236
48 -685854.926 -1576843.870
49 -208261.606 -685854.926
50 -180762.360 -208261.606
51 -528699.952 -180762.360
52 1003307.346 -528699.952
53 -9849.946 1003307.346
54 398402.568 -9849.946
55 1093212.802 398402.568
56 1779267.534 1093212.802
57 1137738.486 1779267.534
58 2191650.424 1137738.486
59 1388478.400 2191650.424
60 NA 1388478.400
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -1703765.203 -1351242.693
[2,] -3600153.537 -1703765.203
[3,] -744792.909 -3600153.537
[4,] 268353.099 -744792.909
[5,] -967670.533 268353.099
[6,] -621709.269 -967670.533
[7,] 280420.255 -621709.269
[8,] -2092034.703 280420.255
[9,] -589810.661 -2092034.703
[10,] -1490518.973 -589810.661
[11,] -2464752.907 -1490518.973
[12,] -410318.313 -2464752.907
[13,] -994530.523 -410318.313
[14,] -1149236.037 -994530.523
[15,] -519295.329 -1149236.037
[16,] -1541743.851 -519295.329
[17,] -400277.263 -1541743.851
[18,] -547553.039 -400277.263
[19,] -304297.785 -547553.039
[20,] -1301421.203 -304297.785
[21,] 126841.209 -1301421.203
[22,] -679947.653 126841.209
[23,] 148188.253 -679947.653
[24,] 326828.687 148188.253
[25,] -207582.783 326828.687
[26,] 1472472.783 -207582.783
[27,] 1155100.501 1472472.783
[28,] 317856.889 1155100.501
[29,] 2330411.147 317856.889
[30,] 1729931.411 2330411.147
[31,] 833561.675 1729931.411
[32,] 3450566.787 833561.675
[33,] 1639176.549 3450566.787
[34,] 1863288.437 1639176.549
[35,] 2504930.123 1863288.437
[36,] 2120587.247 2504930.123
[37,] 3114140.117 2120587.247
[38,] 3457679.150 3114140.117
[39,] 637687.688 3457679.150
[40,] -47773.484 637687.688
[41,] -952613.406 -47773.484
[42,] -959071.672 -952613.406
[43,] -1902896.948 -959071.672
[44,] -1836378.416 -1902896.948
[45,] -2313945.584 -1836378.416
[46,] -1884472.236 -2313945.584
[47,] -1576843.870 -1884472.236
[48,] -685854.926 -1576843.870
[49,] -208261.606 -685854.926
[50,] -180762.360 -208261.606
[51,] -528699.952 -180762.360
[52,] 1003307.346 -528699.952
[53,] -9849.946 1003307.346
[54,] 398402.568 -9849.946
[55,] 1093212.802 398402.568
[56,] 1779267.534 1093212.802
[57,] 1137738.486 1779267.534
[58,] 2191650.424 1137738.486
[59,] 1388478.400 2191650.424
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -1703765.203 -1351242.693
2 -3600153.537 -1703765.203
3 -744792.909 -3600153.537
4 268353.099 -744792.909
5 -967670.533 268353.099
6 -621709.269 -967670.533
7 280420.255 -621709.269
8 -2092034.703 280420.255
9 -589810.661 -2092034.703
10 -1490518.973 -589810.661
11 -2464752.907 -1490518.973
12 -410318.313 -2464752.907
13 -994530.523 -410318.313
14 -1149236.037 -994530.523
15 -519295.329 -1149236.037
16 -1541743.851 -519295.329
17 -400277.263 -1541743.851
18 -547553.039 -400277.263
19 -304297.785 -547553.039
20 -1301421.203 -304297.785
21 126841.209 -1301421.203
22 -679947.653 126841.209
23 148188.253 -679947.653
24 326828.687 148188.253
25 -207582.783 326828.687
26 1472472.783 -207582.783
27 1155100.501 1472472.783
28 317856.889 1155100.501
29 2330411.147 317856.889
30 1729931.411 2330411.147
31 833561.675 1729931.411
32 3450566.787 833561.675
33 1639176.549 3450566.787
34 1863288.437 1639176.549
35 2504930.123 1863288.437
36 2120587.247 2504930.123
37 3114140.117 2120587.247
38 3457679.150 3114140.117
39 637687.688 3457679.150
40 -47773.484 637687.688
41 -952613.406 -47773.484
42 -959071.672 -952613.406
43 -1902896.948 -959071.672
44 -1836378.416 -1902896.948
45 -2313945.584 -1836378.416
46 -1884472.236 -2313945.584
47 -1576843.870 -1884472.236
48 -685854.926 -1576843.870
49 -208261.606 -685854.926
50 -180762.360 -208261.606
51 -528699.952 -180762.360
52 1003307.346 -528699.952
53 -9849.946 1003307.346
54 398402.568 -9849.946
55 1093212.802 398402.568
56 1779267.534 1093212.802
57 1137738.486 1779267.534
58 2191650.424 1137738.486
59 1388478.400 2191650.424
> 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/7ffsq1291458480.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/8ffsq1291458480.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/9qort1291458480.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/10qort1291458480.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/11t7qz1291458480.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/12f7o51291458480.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/13lq3z1291458480.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/14ei2k1291458480.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/15i01q1291458480.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/16dazy1291458480.tab")
+ }
>
> try(system("convert tmp/1jnuh1291458480.ps tmp/1jnuh1291458480.png",intern=TRUE))
character(0)
> try(system("convert tmp/2ufb21291458480.ps tmp/2ufb21291458480.png",intern=TRUE))
character(0)
> try(system("convert tmp/3ufb21291458480.ps tmp/3ufb21291458480.png",intern=TRUE))
character(0)
> try(system("convert tmp/4ufb21291458480.ps tmp/4ufb21291458480.png",intern=TRUE))
character(0)
> try(system("convert tmp/5mosn1291458480.ps tmp/5mosn1291458480.png",intern=TRUE))
character(0)
> try(system("convert tmp/6mosn1291458480.ps tmp/6mosn1291458480.png",intern=TRUE))
character(0)
> try(system("convert tmp/7ffsq1291458480.ps tmp/7ffsq1291458480.png",intern=TRUE))
character(0)
> try(system("convert tmp/8ffsq1291458480.ps tmp/8ffsq1291458480.png",intern=TRUE))
character(0)
> try(system("convert tmp/9qort1291458480.ps tmp/9qort1291458480.png",intern=TRUE))
character(0)
> try(system("convert tmp/10qort1291458480.ps tmp/10qort1291458480.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.791 2.477 4.088