R version 2.9.0 (2009-04-17)
Copyright (C) 2009 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> x <- array(list(921365,0,0,987921,0,0,1132614,0,0,1332224,0,0,1418133,0,0,1411549,0,0,1695920,0,0,1636173,0,0,1539653,0,0,1395314,0,0,1127575,0,0,1036076,0,0,989236,0,0,1008380,0,0,1207763,0,0,1368839,0,0,1469798,0,0,1498721,0,0,1761769,0,0,1653214,0,0,1599104,0,0,1421179,0,0,1163995,0,0,1037735,0,0,1015407,0,0,1039210,0,0,1258049,0,0,1469445,0,0,1552346,0,0,1549144,0,0,1785895,0,0,1662335,0,0,1629440,0,0,1467430,0,0,1202209,0,0,1076982,0,0,1039367,1,0,1063449,1,0,1335135,1,0,1491602,1,0,1591972,1,0,1641248,1,0,1898849,1,0,1798580,1,0,1762444,1,0,1622044,1,0,1368955,1,0,1262973,1,0,1195650,1,0,1269530,1,0,1479279,1,0,1607819,1,0,1712466,1,0,1721766,1,0,1949843,1,0,1821326,1,0,1757802,1,1,1590367,1,1,1260647,1,1,1149235,1,1),dim=c(3,60),dimnames=list(c('Y','X','Crisis'),1:60))
> y <- array(NA,dim=c(3,60),dimnames=list(c('Y','X','Crisis'),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 = '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 Crisis M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 921365 0 0 1 0 0 0 0 0 0 0 0 0 0 1
2 987921 0 0 0 1 0 0 0 0 0 0 0 0 0 2
3 1132614 0 0 0 0 1 0 0 0 0 0 0 0 0 3
4 1332224 0 0 0 0 0 1 0 0 0 0 0 0 0 4
5 1418133 0 0 0 0 0 0 1 0 0 0 0 0 0 5
6 1411549 0 0 0 0 0 0 0 1 0 0 0 0 0 6
7 1695920 0 0 0 0 0 0 0 0 1 0 0 0 0 7
8 1636173 0 0 0 0 0 0 0 0 0 1 0 0 0 8
9 1539653 0 0 0 0 0 0 0 0 0 0 1 0 0 9
10 1395314 0 0 0 0 0 0 0 0 0 0 0 1 0 10
11 1127575 0 0 0 0 0 0 0 0 0 0 0 0 1 11
12 1036076 0 0 0 0 0 0 0 0 0 0 0 0 0 12
13 989236 0 0 1 0 0 0 0 0 0 0 0 0 0 13
14 1008380 0 0 0 1 0 0 0 0 0 0 0 0 0 14
15 1207763 0 0 0 0 1 0 0 0 0 0 0 0 0 15
16 1368839 0 0 0 0 0 1 0 0 0 0 0 0 0 16
17 1469798 0 0 0 0 0 0 1 0 0 0 0 0 0 17
18 1498721 0 0 0 0 0 0 0 1 0 0 0 0 0 18
19 1761769 0 0 0 0 0 0 0 0 1 0 0 0 0 19
20 1653214 0 0 0 0 0 0 0 0 0 1 0 0 0 20
21 1599104 0 0 0 0 0 0 0 0 0 0 1 0 0 21
22 1421179 0 0 0 0 0 0 0 0 0 0 0 1 0 22
23 1163995 0 0 0 0 0 0 0 0 0 0 0 0 1 23
24 1037735 0 0 0 0 0 0 0 0 0 0 0 0 0 24
25 1015407 0 0 1 0 0 0 0 0 0 0 0 0 0 25
26 1039210 0 0 0 1 0 0 0 0 0 0 0 0 0 26
27 1258049 0 0 0 0 1 0 0 0 0 0 0 0 0 27
28 1469445 0 0 0 0 0 1 0 0 0 0 0 0 0 28
29 1552346 0 0 0 0 0 0 1 0 0 0 0 0 0 29
30 1549144 0 0 0 0 0 0 0 1 0 0 0 0 0 30
31 1785895 0 0 0 0 0 0 0 0 1 0 0 0 0 31
32 1662335 0 0 0 0 0 0 0 0 0 1 0 0 0 32
33 1629440 0 0 0 0 0 0 0 0 0 0 1 0 0 33
34 1467430 0 0 0 0 0 0 0 0 0 0 0 1 0 34
35 1202209 0 0 0 0 0 0 0 0 0 0 0 0 1 35
36 1076982 0 0 0 0 0 0 0 0 0 0 0 0 0 36
37 1039367 1 0 1 0 0 0 0 0 0 0 0 0 0 37
38 1063449 1 0 0 1 0 0 0 0 0 0 0 0 0 38
39 1335135 1 0 0 0 1 0 0 0 0 0 0 0 0 39
40 1491602 1 0 0 0 0 1 0 0 0 0 0 0 0 40
41 1591972 1 0 0 0 0 0 1 0 0 0 0 0 0 41
42 1641248 1 0 0 0 0 0 0 1 0 0 0 0 0 42
43 1898849 1 0 0 0 0 0 0 0 1 0 0 0 0 43
44 1798580 1 0 0 0 0 0 0 0 0 1 0 0 0 44
45 1762444 1 0 0 0 0 0 0 0 0 0 1 0 0 45
46 1622044 1 0 0 0 0 0 0 0 0 0 0 1 0 46
47 1368955 1 0 0 0 0 0 0 0 0 0 0 0 1 47
48 1262973 1 0 0 0 0 0 0 0 0 0 0 0 0 48
49 1195650 1 0 1 0 0 0 0 0 0 0 0 0 0 49
50 1269530 1 0 0 1 0 0 0 0 0 0 0 0 0 50
51 1479279 1 0 0 0 1 0 0 0 0 0 0 0 0 51
52 1607819 1 0 0 0 0 1 0 0 0 0 0 0 0 52
53 1712466 1 0 0 0 0 0 1 0 0 0 0 0 0 53
54 1721766 1 0 0 0 0 0 0 1 0 0 0 0 0 54
55 1949843 1 0 0 0 0 0 0 0 1 0 0 0 0 55
56 1821326 1 0 0 0 0 0 0 0 0 1 0 0 0 56
57 1757802 1 1 0 0 0 0 0 0 0 0 1 0 0 57
58 1590367 1 1 0 0 0 0 0 0 0 0 0 1 0 58
59 1260647 1 1 0 0 0 0 0 0 0 0 0 0 1 59
60 1149235 1 1 0 0 0 0 0 0 0 0 0 0 0 60
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X Crisis M1 M2 M3
952334 52433 -87551 -49993 -12856 191658
M4 M5 M6 M7 M8 M9
358721 449322 460509 710123 601638 558155
M10 M11 t
395378 116432 4356
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-93976.5 -19940.0 -518.6 21348.7 60715.7
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 952333.9 20235.1 47.063 < 2e-16 ***
X 52433.0 18171.7 2.885 0.005982 **
Crisis -87550.6 21113.6 -4.147 0.000147 ***
M1 -49993.2 22746.4 -2.198 0.033146 *
M2 -12855.9 22634.5 -0.568 0.572874
M3 191658.5 22534.9 8.505 6.42e-11 ***
M4 358720.6 22447.7 15.980 < 2e-16 ***
M5 449322.2 22373.1 20.083 < 2e-16 ***
M6 460509.2 22311.2 20.640 < 2e-16 ***
M7 710123.1 22262.1 31.898 < 2e-16 ***
M8 601637.9 22226.0 27.069 < 2e-16 ***
M9 558155.3 21865.5 25.527 < 2e-16 ***
M10 395377.9 21832.5 18.110 < 2e-16 ***
M11 116431.6 21812.6 5.338 2.96e-06 ***
t 4355.6 537.5 8.103 2.44e-10 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 34480 on 45 degrees of freedom
Multiple R-squared: 0.9878, Adjusted R-squared: 0.984
F-statistic: 260.5 on 14 and 45 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.219946285 0.439892570 0.7800537
[2,] 0.113648088 0.227296175 0.8863519
[3,] 0.122583450 0.245166899 0.8774166
[4,] 0.062040501 0.124081002 0.9379595
[5,] 0.040159816 0.080319633 0.9598402
[6,] 0.021111595 0.042223190 0.9788884
[7,] 0.027310798 0.054621597 0.9726892
[8,] 0.015269047 0.030538093 0.9847310
[9,] 0.010153440 0.020306879 0.9898466
[10,] 0.006382392 0.012764784 0.9936176
[11,] 0.020006594 0.040013188 0.9799934
[12,] 0.028362842 0.056725683 0.9716372
[13,] 0.022156991 0.044313982 0.9778430
[14,] 0.020094151 0.040188302 0.9799058
[15,] 0.051807484 0.103614968 0.9481925
[16,] 0.030169966 0.060339932 0.9698300
[17,] 0.017205114 0.034410228 0.9827949
[18,] 0.009283343 0.018566686 0.9907167
[19,] 0.005689161 0.011378323 0.9943108
[20,] 0.003679164 0.007358329 0.9963208
[21,] 0.018546918 0.037093837 0.9814531
[22,] 0.059408493 0.118816986 0.9405915
[23,] 0.049096678 0.098193355 0.9509033
[24,] 0.068228624 0.136457248 0.9317714
[25,] 0.073578493 0.147156985 0.9264215
> postscript(file="/var/www/html/rcomp/tmp/1ueon1261311068.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> plot(x[,1], type='l', main='Actuals and Interpolation', ylab='value of Actuals and Interpolation (dots)', xlab='time or index')
> points(x[,1]-mysum$resid)
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/23uuj1261311068.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> plot(mysum$resid, type='b', pch=19, main='Residuals', ylab='value of Residuals', xlab='time or index')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/3gxsz1261311068.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> hist(mysum$resid, main='Residual Histogram', xlab='values of Residuals')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/4m3or1261311068.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> densityplot(~mysum$resid,col='black',main='Residual Density Plot', xlab='values of Residuals')
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/5ejvw1261311068.ps",horizontal=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
14668.6656 39731.6656 -24445.3344 3746.8656 -5301.3344 -27427.9344
7 8 9 10 11 12
2973.4656 47356.0656 -10037.0494 4045.7506 10897.3506 31474.3506
13 14 15 16 17 18
30271.9272 7922.9272 -1564.0728 -11905.8728 -5904.0728 7476.3272
19 20 21 22 23 24
16554.7272 12129.3272 -2853.7878 -22356.9878 -4950.3878 -19134.3878
25 26 27 28 29 30
4175.1889 -13514.8111 -3545.8111 36432.3889 24376.1889 5631.5889
31 32 33 34 35 36
-11587.0111 -31017.4111 -24785.5261 -28373.7261 -19004.1261 -32155.1261
37 38 39 40 41 42
-76565.5217 -93976.5217 -31160.5217 -46111.3217 -40698.5217 -6965.1217
43 44 45 46 47 48
-3333.7217 526.8783 3517.7633 21539.5633 43041.1633 49135.1633
49 50 51 52 53 54
27449.7400 59836.7400 60715.7400 17837.9400 27527.7400 21285.1400
55 56 57 58 59 60
-4607.4600 -28994.8600 34158.6000 25145.4000 -29984.0000 -29320.0000
> postscript(file="/var/www/html/rcomp/tmp/6acis1261311068.ps",horizontal=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 14668.6656 NA
1 39731.6656 14668.6656
2 -24445.3344 39731.6656
3 3746.8656 -24445.3344
4 -5301.3344 3746.8656
5 -27427.9344 -5301.3344
6 2973.4656 -27427.9344
7 47356.0656 2973.4656
8 -10037.0494 47356.0656
9 4045.7506 -10037.0494
10 10897.3506 4045.7506
11 31474.3506 10897.3506
12 30271.9272 31474.3506
13 7922.9272 30271.9272
14 -1564.0728 7922.9272
15 -11905.8728 -1564.0728
16 -5904.0728 -11905.8728
17 7476.3272 -5904.0728
18 16554.7272 7476.3272
19 12129.3272 16554.7272
20 -2853.7878 12129.3272
21 -22356.9878 -2853.7878
22 -4950.3878 -22356.9878
23 -19134.3878 -4950.3878
24 4175.1889 -19134.3878
25 -13514.8111 4175.1889
26 -3545.8111 -13514.8111
27 36432.3889 -3545.8111
28 24376.1889 36432.3889
29 5631.5889 24376.1889
30 -11587.0111 5631.5889
31 -31017.4111 -11587.0111
32 -24785.5261 -31017.4111
33 -28373.7261 -24785.5261
34 -19004.1261 -28373.7261
35 -32155.1261 -19004.1261
36 -76565.5217 -32155.1261
37 -93976.5217 -76565.5217
38 -31160.5217 -93976.5217
39 -46111.3217 -31160.5217
40 -40698.5217 -46111.3217
41 -6965.1217 -40698.5217
42 -3333.7217 -6965.1217
43 526.8783 -3333.7217
44 3517.7633 526.8783
45 21539.5633 3517.7633
46 43041.1633 21539.5633
47 49135.1633 43041.1633
48 27449.7400 49135.1633
49 59836.7400 27449.7400
50 60715.7400 59836.7400
51 17837.9400 60715.7400
52 27527.7400 17837.9400
53 21285.1400 27527.7400
54 -4607.4600 21285.1400
55 -28994.8600 -4607.4600
56 34158.6000 -28994.8600
57 25145.4000 34158.6000
58 -29984.0000 25145.4000
59 -29320.0000 -29984.0000
60 NA -29320.0000
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 39731.6656 14668.6656
[2,] -24445.3344 39731.6656
[3,] 3746.8656 -24445.3344
[4,] -5301.3344 3746.8656
[5,] -27427.9344 -5301.3344
[6,] 2973.4656 -27427.9344
[7,] 47356.0656 2973.4656
[8,] -10037.0494 47356.0656
[9,] 4045.7506 -10037.0494
[10,] 10897.3506 4045.7506
[11,] 31474.3506 10897.3506
[12,] 30271.9272 31474.3506
[13,] 7922.9272 30271.9272
[14,] -1564.0728 7922.9272
[15,] -11905.8728 -1564.0728
[16,] -5904.0728 -11905.8728
[17,] 7476.3272 -5904.0728
[18,] 16554.7272 7476.3272
[19,] 12129.3272 16554.7272
[20,] -2853.7878 12129.3272
[21,] -22356.9878 -2853.7878
[22,] -4950.3878 -22356.9878
[23,] -19134.3878 -4950.3878
[24,] 4175.1889 -19134.3878
[25,] -13514.8111 4175.1889
[26,] -3545.8111 -13514.8111
[27,] 36432.3889 -3545.8111
[28,] 24376.1889 36432.3889
[29,] 5631.5889 24376.1889
[30,] -11587.0111 5631.5889
[31,] -31017.4111 -11587.0111
[32,] -24785.5261 -31017.4111
[33,] -28373.7261 -24785.5261
[34,] -19004.1261 -28373.7261
[35,] -32155.1261 -19004.1261
[36,] -76565.5217 -32155.1261
[37,] -93976.5217 -76565.5217
[38,] -31160.5217 -93976.5217
[39,] -46111.3217 -31160.5217
[40,] -40698.5217 -46111.3217
[41,] -6965.1217 -40698.5217
[42,] -3333.7217 -6965.1217
[43,] 526.8783 -3333.7217
[44,] 3517.7633 526.8783
[45,] 21539.5633 3517.7633
[46,] 43041.1633 21539.5633
[47,] 49135.1633 43041.1633
[48,] 27449.7400 49135.1633
[49,] 59836.7400 27449.7400
[50,] 60715.7400 59836.7400
[51,] 17837.9400 60715.7400
[52,] 27527.7400 17837.9400
[53,] 21285.1400 27527.7400
[54,] -4607.4600 21285.1400
[55,] -28994.8600 -4607.4600
[56,] 34158.6000 -28994.8600
[57,] 25145.4000 34158.6000
[58,] -29984.0000 25145.4000
[59,] -29320.0000 -29984.0000
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 39731.6656 14668.6656
2 -24445.3344 39731.6656
3 3746.8656 -24445.3344
4 -5301.3344 3746.8656
5 -27427.9344 -5301.3344
6 2973.4656 -27427.9344
7 47356.0656 2973.4656
8 -10037.0494 47356.0656
9 4045.7506 -10037.0494
10 10897.3506 4045.7506
11 31474.3506 10897.3506
12 30271.9272 31474.3506
13 7922.9272 30271.9272
14 -1564.0728 7922.9272
15 -11905.8728 -1564.0728
16 -5904.0728 -11905.8728
17 7476.3272 -5904.0728
18 16554.7272 7476.3272
19 12129.3272 16554.7272
20 -2853.7878 12129.3272
21 -22356.9878 -2853.7878
22 -4950.3878 -22356.9878
23 -19134.3878 -4950.3878
24 4175.1889 -19134.3878
25 -13514.8111 4175.1889
26 -3545.8111 -13514.8111
27 36432.3889 -3545.8111
28 24376.1889 36432.3889
29 5631.5889 24376.1889
30 -11587.0111 5631.5889
31 -31017.4111 -11587.0111
32 -24785.5261 -31017.4111
33 -28373.7261 -24785.5261
34 -19004.1261 -28373.7261
35 -32155.1261 -19004.1261
36 -76565.5217 -32155.1261
37 -93976.5217 -76565.5217
38 -31160.5217 -93976.5217
39 -46111.3217 -31160.5217
40 -40698.5217 -46111.3217
41 -6965.1217 -40698.5217
42 -3333.7217 -6965.1217
43 526.8783 -3333.7217
44 3517.7633 526.8783
45 21539.5633 3517.7633
46 43041.1633 21539.5633
47 49135.1633 43041.1633
48 27449.7400 49135.1633
49 59836.7400 27449.7400
50 60715.7400 59836.7400
51 17837.9400 60715.7400
52 27527.7400 17837.9400
53 21285.1400 27527.7400
54 -4607.4600 21285.1400
55 -28994.8600 -4607.4600
56 34158.6000 -28994.8600
57 25145.4000 34158.6000
58 -29984.0000 25145.4000
59 -29320.0000 -29984.0000
> plot(z,main=paste('Residual Lag plot, lowess, and regression line'), ylab='values of Residuals', xlab='lagged values of Residuals')
> lines(lowess(z))
> abline(lm(z))
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/7xwri1261311068.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> acf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Autocorrelation Function')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/8lyor1261311068.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> pacf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Partial Autocorrelation Function')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/92c3w1261311068.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> opar <- par(mfrow = c(2,2), oma = c(0, 0, 1.1, 0))
> plot(mylm, las = 1, sub='Residual Diagnostics')
> par(opar)
> dev.off()
null device
1
> if (n > n25) {
+ postscript(file="/var/www/html/rcomp/tmp/10ij031261311068.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
+ plot(kp3:nmkm3,gqarr[,2], main='Goldfeld-Quandt test',ylab='2-sided p-value',xlab='breakpoint')
+ grid()
+ dev.off()
+ }
null device
1
>
> #Note: the /var/www/html/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/html/rcomp/createtable")
>
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Estimated Regression Equation', 1, TRUE)
> a<-table.row.end(a)
> myeq <- colnames(x)[1]
> myeq <- paste(myeq, '[t] = ', sep='')
> for (i in 1:k){
+ if (mysum$coefficients[i,1] > 0) myeq <- paste(myeq, '+', '')
+ myeq <- paste(myeq, mysum$coefficients[i,1], sep=' ')
+ if (rownames(mysum$coefficients)[i] != '(Intercept)') {
+ myeq <- paste(myeq, rownames(mysum$coefficients)[i], sep='')
+ if (rownames(mysum$coefficients)[i] != 't') myeq <- paste(myeq, '[t]', sep='')
+ }
+ }
> myeq <- paste(myeq, ' + e[t]')
> a<-table.row.start(a)
> a<-table.element(a, myeq)
> a<-table.row.end(a)
> a<-table.end(a)
> table.save(a,file="/var/www/html/rcomp/tmp/117ydk1261311068.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a,hyperlink('http://www.xycoon.com/ols1.htm','Multiple Linear Regression - Ordinary Least Squares',''), 6, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a,'Variable',header=TRUE)
> a<-table.element(a,'Parameter',header=TRUE)
> a<-table.element(a,'S.D.',header=TRUE)
> a<-table.element(a,'T-STAT
H0: parameter = 0',header=TRUE)
> a<-table.element(a,'2-tail p-value',header=TRUE)
> a<-table.element(a,'1-tail p-value',header=TRUE)
> a<-table.row.end(a)
> for (i in 1:k){
+ a<-table.row.start(a)
+ a<-table.element(a,rownames(mysum$coefficients)[i],header=TRUE)
+ a<-table.element(a,mysum$coefficients[i,1])
+ a<-table.element(a, round(mysum$coefficients[i,2],6))
+ a<-table.element(a, round(mysum$coefficients[i,3],4))
+ a<-table.element(a, round(mysum$coefficients[i,4],6))
+ a<-table.element(a, round(mysum$coefficients[i,4]/2,6))
+ a<-table.row.end(a)
+ }
> a<-table.end(a)
> table.save(a,file="/var/www/html/rcomp/tmp/12tc6y1261311068.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Regression Statistics', 2, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple R',1,TRUE)
> a<-table.element(a, sqrt(mysum$r.squared))
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'R-squared',1,TRUE)
> a<-table.element(a, mysum$r.squared)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Adjusted R-squared',1,TRUE)
> a<-table.element(a, mysum$adj.r.squared)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (value)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[1])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (DF numerator)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[2])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (DF denominator)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[3])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'p-value',1,TRUE)
> a<-table.element(a, 1-pf(mysum$fstatistic[1],mysum$fstatistic[2],mysum$fstatistic[3]))
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Residual Statistics', 2, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Residual Standard Deviation',1,TRUE)
> a<-table.element(a, mysum$sigma)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Sum Squared Residuals',1,TRUE)
> a<-table.element(a, sum(myerror*myerror))
> a<-table.row.end(a)
> a<-table.end(a)
> table.save(a,file="/var/www/html/rcomp/tmp/13c4nu1261311068.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Actuals, Interpolation, and Residuals', 4, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Time or Index', 1, TRUE)
> a<-table.element(a, 'Actuals', 1, TRUE)
> a<-table.element(a, 'Interpolation
Forecast', 1, TRUE)
> a<-table.element(a, 'Residuals
Prediction Error', 1, TRUE)
> a<-table.row.end(a)
> for (i in 1:n) {
+ a<-table.row.start(a)
+ a<-table.element(a,i, 1, TRUE)
+ a<-table.element(a,x[i])
+ a<-table.element(a,x[i]-mysum$resid[i])
+ a<-table.element(a,mysum$resid[i])
+ a<-table.row.end(a)
+ }
> a<-table.end(a)
> table.save(a,file="/var/www/html/rcomp/tmp/148row1261311068.tab")
> if (n > n25) {
+ a<-table.start()
+ a<-table.row.start(a)
+ a<-table.element(a,'Goldfeld-Quandt test for Heteroskedasticity',4,TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'p-values',header=TRUE)
+ a<-table.element(a,'Alternative Hypothesis',3,header=TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'breakpoint index',header=TRUE)
+ a<-table.element(a,'greater',header=TRUE)
+ a<-table.element(a,'2-sided',header=TRUE)
+ a<-table.element(a,'less',header=TRUE)
+ a<-table.row.end(a)
+ for (mypoint in kp3:nmkm3) {
+ a<-table.row.start(a)
+ a<-table.element(a,mypoint,header=TRUE)
+ a<-table.element(a,gqarr[mypoint-kp3+1,1])
+ a<-table.element(a,gqarr[mypoint-kp3+1,2])
+ a<-table.element(a,gqarr[mypoint-kp3+1,3])
+ a<-table.row.end(a)
+ }
+ a<-table.end(a)
+ table.save(a,file="/var/www/html/rcomp/tmp/15gubv1261311068.tab")
+ a<-table.start()
+ a<-table.row.start(a)
+ a<-table.element(a,'Meta Analysis of Goldfeld-Quandt test for Heteroskedasticity',4,TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'Description',header=TRUE)
+ a<-table.element(a,'# significant tests',header=TRUE)
+ a<-table.element(a,'% significant tests',header=TRUE)
+ a<-table.element(a,'OK/NOK',header=TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'1% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant1)
+ a<-table.element(a,numsignificant1/numgqtests)
+ if (numsignificant1/numgqtests < 0.01) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'5% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant5)
+ a<-table.element(a,numsignificant5/numgqtests)
+ if (numsignificant5/numgqtests < 0.05) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'10% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant10)
+ a<-table.element(a,numsignificant10/numgqtests)
+ if (numsignificant10/numgqtests < 0.1) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.end(a)
+ table.save(a,file="/var/www/html/rcomp/tmp/16gp501261311068.tab")
+ }
>
> try(system("convert tmp/1ueon1261311068.ps tmp/1ueon1261311068.png",intern=TRUE))
character(0)
> try(system("convert tmp/23uuj1261311068.ps tmp/23uuj1261311068.png",intern=TRUE))
character(0)
> try(system("convert tmp/3gxsz1261311068.ps tmp/3gxsz1261311068.png",intern=TRUE))
character(0)
> try(system("convert tmp/4m3or1261311068.ps tmp/4m3or1261311068.png",intern=TRUE))
character(0)
> try(system("convert tmp/5ejvw1261311068.ps tmp/5ejvw1261311068.png",intern=TRUE))
character(0)
> try(system("convert tmp/6acis1261311068.ps tmp/6acis1261311068.png",intern=TRUE))
character(0)
> try(system("convert tmp/7xwri1261311068.ps tmp/7xwri1261311068.png",intern=TRUE))
character(0)
> try(system("convert tmp/8lyor1261311068.ps tmp/8lyor1261311068.png",intern=TRUE))
character(0)
> try(system("convert tmp/92c3w1261311068.ps tmp/92c3w1261311068.png",intern=TRUE))
character(0)
> try(system("convert tmp/10ij031261311068.ps tmp/10ij031261311068.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
2.389 1.543 4.800