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(286602,0,283042,0,276687,0,277915,0,277128,0,277103,0,275037,0,270150,0,267140,0,264993,0,287259,0,291186,0,292300,0,288186,0,281477,0,282656,0,280190,0,280408,0,276836,0,275216,0,274352,0,271311,0,289802,0,290726,0,292300,0,278506,0,269826,0,265861,0,269034,0,264176,0,255198,0,253353,0,246057,0,235372,0,258556,0,260993,0,254663,0,250643,0,243422,0,247105,0,248541,0,245039,0,237080,0,237085,0,225554,0,226839,1,247934,1,248333,1,246969,1,245098,1,246263,1,255765,1,264319,1,268347,1,273046,1,273963,1,267430,1,271993,1,292710,1,295881,1),dim=c(2,60),dimnames=list(c('nwwmb','dummy_variable'),1:60))
> y <- array(NA,dim=c(2,60),dimnames=list(c('nwwmb','dummy_variable'),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
nwwmb dummy_variable M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 286602 0 1 0 0 0 0 0 0 0 0 0 0
2 283042 0 0 1 0 0 0 0 0 0 0 0 0
3 276687 0 0 0 1 0 0 0 0 0 0 0 0
4 277915 0 0 0 0 1 0 0 0 0 0 0 0
5 277128 0 0 0 0 0 1 0 0 0 0 0 0
6 277103 0 0 0 0 0 0 1 0 0 0 0 0
7 275037 0 0 0 0 0 0 0 1 0 0 0 0
8 270150 0 0 0 0 0 0 0 0 1 0 0 0
9 267140 0 0 0 0 0 0 0 0 0 1 0 0
10 264993 0 0 0 0 0 0 0 0 0 0 1 0
11 287259 0 0 0 0 0 0 0 0 0 0 0 1
12 291186 0 0 0 0 0 0 0 0 0 0 0 0
13 292300 0 1 0 0 0 0 0 0 0 0 0 0
14 288186 0 0 1 0 0 0 0 0 0 0 0 0
15 281477 0 0 0 1 0 0 0 0 0 0 0 0
16 282656 0 0 0 0 1 0 0 0 0 0 0 0
17 280190 0 0 0 0 0 1 0 0 0 0 0 0
18 280408 0 0 0 0 0 0 1 0 0 0 0 0
19 276836 0 0 0 0 0 0 0 1 0 0 0 0
20 275216 0 0 0 0 0 0 0 0 1 0 0 0
21 274352 0 0 0 0 0 0 0 0 0 1 0 0
22 271311 0 0 0 0 0 0 0 0 0 0 1 0
23 289802 0 0 0 0 0 0 0 0 0 0 0 1
24 290726 0 0 0 0 0 0 0 0 0 0 0 0
25 292300 0 1 0 0 0 0 0 0 0 0 0 0
26 278506 0 0 1 0 0 0 0 0 0 0 0 0
27 269826 0 0 0 1 0 0 0 0 0 0 0 0
28 265861 0 0 0 0 1 0 0 0 0 0 0 0
29 269034 0 0 0 0 0 1 0 0 0 0 0 0
30 264176 0 0 0 0 0 0 1 0 0 0 0 0
31 255198 0 0 0 0 0 0 0 1 0 0 0 0
32 253353 0 0 0 0 0 0 0 0 1 0 0 0
33 246057 0 0 0 0 0 0 0 0 0 1 0 0
34 235372 0 0 0 0 0 0 0 0 0 0 1 0
35 258556 0 0 0 0 0 0 0 0 0 0 0 1
36 260993 0 0 0 0 0 0 0 0 0 0 0 0
37 254663 0 1 0 0 0 0 0 0 0 0 0 0
38 250643 0 0 1 0 0 0 0 0 0 0 0 0
39 243422 0 0 0 1 0 0 0 0 0 0 0 0
40 247105 0 0 0 0 1 0 0 0 0 0 0 0
41 248541 0 0 0 0 0 1 0 0 0 0 0 0
42 245039 0 0 0 0 0 0 1 0 0 0 0 0
43 237080 0 0 0 0 0 0 0 1 0 0 0 0
44 237085 0 0 0 0 0 0 0 0 1 0 0 0
45 225554 0 0 0 0 0 0 0 0 0 1 0 0
46 226839 1 0 0 0 0 0 0 0 0 0 1 0
47 247934 1 0 0 0 0 0 0 0 0 0 0 1
48 248333 1 0 0 0 0 0 0 0 0 0 0 0
49 246969 1 1 0 0 0 0 0 0 0 0 0 0
50 245098 1 0 1 0 0 0 0 0 0 0 0 0
51 246263 1 0 0 1 0 0 0 0 0 0 0 0
52 255765 1 0 0 0 1 0 0 0 0 0 0 0
53 264319 1 0 0 0 0 1 0 0 0 0 0 0
54 268347 1 0 0 0 0 0 1 0 0 0 0 0
55 273046 1 0 0 0 0 0 0 1 0 0 0 0
56 273963 1 0 0 0 0 0 0 0 1 0 0 0
57 267430 1 0 0 0 0 0 0 0 0 1 0 0
58 271993 1 0 0 0 0 0 0 0 0 0 1 0
59 292710 1 0 0 0 0 0 0 0 0 0 0 1
60 295881 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) dummy_variable M1 M2 M3
280316 -7230 -4303 -9775 -15335
M4 M5 M6 M7 M8
-13009 -11027 -11855 -15430 -16916
M9 M10 M11
-22763 -23322 -2172
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-31999 -19389 7902 11948 22795
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 280316 8458 33.142 <2e-16 ***
dummy_variable -7230 5553 -1.302 0.1993
M1 -4303 11595 -0.371 0.7122
M2 -9775 11595 -0.843 0.4035
M3 -15335 11595 -1.323 0.1924
M4 -13009 11595 -1.122 0.2676
M5 -11027 11595 -0.951 0.3465
M6 -11855 11595 -1.022 0.3118
M7 -15430 11595 -1.331 0.1897
M8 -16916 11595 -1.459 0.1512
M9 -22763 11595 -1.963 0.0556 .
M10 -23322 11542 -2.021 0.0490 *
M11 -2172 11542 -0.188 0.8516
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 18250 on 47 degrees of freedom
Multiple R-squared: 0.1792, Adjusted R-squared: -0.03038
F-statistic: 0.855 on 12 and 47 DF, p-value: 0.5958
> 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,] 1.320034e-02 2.640068e-02 0.9867997
[2,] 2.762821e-03 5.525642e-03 0.9972372
[3,] 5.905518e-04 1.181104e-03 0.9994094
[4,] 1.070981e-04 2.141963e-04 0.9998929
[5,] 3.194514e-05 6.389028e-05 0.9999681
[6,] 1.850743e-05 3.701486e-05 0.9999815
[7,] 8.629697e-06 1.725939e-05 0.9999914
[8,] 2.226819e-06 4.453638e-06 0.9999978
[9,] 5.212992e-07 1.042598e-06 0.9999995
[10,] 3.114686e-07 6.229372e-07 0.9999997
[11,] 4.894650e-07 9.789300e-07 0.9999995
[12,] 1.335833e-06 2.671666e-06 0.9999987
[13,] 1.112481e-05 2.224961e-05 0.9999889
[14,] 1.274867e-05 2.549734e-05 0.9999873
[15,] 3.481971e-05 6.963942e-05 0.9999652
[16,] 2.077859e-04 4.155717e-04 0.9997922
[17,] 4.529762e-04 9.059524e-04 0.9995470
[18,] 1.502465e-03 3.004930e-03 0.9984975
[19,] 6.990953e-03 1.398191e-02 0.9930090
[20,] 1.295849e-02 2.591697e-02 0.9870415
[21,] 1.849995e-02 3.699990e-02 0.9815000
[22,] 4.523922e-02 9.047844e-02 0.9547608
[23,] 7.409043e-02 1.481809e-01 0.9259096
[24,] 9.411168e-02 1.882234e-01 0.9058883
[25,] 9.439153e-02 1.887831e-01 0.9056085
[26,] 8.149561e-02 1.629912e-01 0.9185044
[27,] 6.414406e-02 1.282881e-01 0.9358559
[28,] 4.441487e-02 8.882975e-02 0.9555851
[29,] 2.528704e-02 5.057409e-02 0.9747130
> postscript(file="/var/www/html/rcomp/tmp/1hje01258742239.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/2l78t1258742239.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/3vo3s1258742239.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/4j86l1258742239.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/5u8co1258742239.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
10589.2963 12501.0963 11706.0963 10608.6963 7839.6963 8642.4963
7 8 9 10 11 12
10151.6963 6750.6963 9587.4963 7999.5926 9114.9926 10870.3926
13 14 15 16 17 18
16287.2963 17645.0963 16496.0963 15349.6963 10901.6963 11947.4963
19 20 21 22 23 24
11950.6963 11816.6963 16799.4963 14317.5926 11657.9926 10410.3926
25 26 27 28 29 30
16287.2963 7965.0963 4845.0963 -1445.3037 -254.3037 -4284.5037
31 32 33 34 35 36
-9687.3037 -10046.3037 -11495.5037 -21621.4074 -19588.0074 -19322.6074
37 38 39 40 41 42
-21349.7037 -19897.9037 -21558.9037 -20201.3037 -20747.3037 -23421.5037
43 44 45 46 47 48
-27805.3037 -26314.3037 -31998.5037 -22924.8889 -22980.4889 -24753.0889
49 50 51 52 53 54
-21814.1852 -18213.3852 -11488.3852 -4311.7852 2260.2148 7116.0148
55 56 57 58 59 60
15390.2148 17793.2148 17107.0148 22229.1111 21795.5111 22794.9111
> postscript(file="/var/www/html/rcomp/tmp/6vwoe1258742239.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 10589.2963 NA
1 12501.0963 10589.2963
2 11706.0963 12501.0963
3 10608.6963 11706.0963
4 7839.6963 10608.6963
5 8642.4963 7839.6963
6 10151.6963 8642.4963
7 6750.6963 10151.6963
8 9587.4963 6750.6963
9 7999.5926 9587.4963
10 9114.9926 7999.5926
11 10870.3926 9114.9926
12 16287.2963 10870.3926
13 17645.0963 16287.2963
14 16496.0963 17645.0963
15 15349.6963 16496.0963
16 10901.6963 15349.6963
17 11947.4963 10901.6963
18 11950.6963 11947.4963
19 11816.6963 11950.6963
20 16799.4963 11816.6963
21 14317.5926 16799.4963
22 11657.9926 14317.5926
23 10410.3926 11657.9926
24 16287.2963 10410.3926
25 7965.0963 16287.2963
26 4845.0963 7965.0963
27 -1445.3037 4845.0963
28 -254.3037 -1445.3037
29 -4284.5037 -254.3037
30 -9687.3037 -4284.5037
31 -10046.3037 -9687.3037
32 -11495.5037 -10046.3037
33 -21621.4074 -11495.5037
34 -19588.0074 -21621.4074
35 -19322.6074 -19588.0074
36 -21349.7037 -19322.6074
37 -19897.9037 -21349.7037
38 -21558.9037 -19897.9037
39 -20201.3037 -21558.9037
40 -20747.3037 -20201.3037
41 -23421.5037 -20747.3037
42 -27805.3037 -23421.5037
43 -26314.3037 -27805.3037
44 -31998.5037 -26314.3037
45 -22924.8889 -31998.5037
46 -22980.4889 -22924.8889
47 -24753.0889 -22980.4889
48 -21814.1852 -24753.0889
49 -18213.3852 -21814.1852
50 -11488.3852 -18213.3852
51 -4311.7852 -11488.3852
52 2260.2148 -4311.7852
53 7116.0148 2260.2148
54 15390.2148 7116.0148
55 17793.2148 15390.2148
56 17107.0148 17793.2148
57 22229.1111 17107.0148
58 21795.5111 22229.1111
59 22794.9111 21795.5111
60 NA 22794.9111
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 12501.0963 10589.2963
[2,] 11706.0963 12501.0963
[3,] 10608.6963 11706.0963
[4,] 7839.6963 10608.6963
[5,] 8642.4963 7839.6963
[6,] 10151.6963 8642.4963
[7,] 6750.6963 10151.6963
[8,] 9587.4963 6750.6963
[9,] 7999.5926 9587.4963
[10,] 9114.9926 7999.5926
[11,] 10870.3926 9114.9926
[12,] 16287.2963 10870.3926
[13,] 17645.0963 16287.2963
[14,] 16496.0963 17645.0963
[15,] 15349.6963 16496.0963
[16,] 10901.6963 15349.6963
[17,] 11947.4963 10901.6963
[18,] 11950.6963 11947.4963
[19,] 11816.6963 11950.6963
[20,] 16799.4963 11816.6963
[21,] 14317.5926 16799.4963
[22,] 11657.9926 14317.5926
[23,] 10410.3926 11657.9926
[24,] 16287.2963 10410.3926
[25,] 7965.0963 16287.2963
[26,] 4845.0963 7965.0963
[27,] -1445.3037 4845.0963
[28,] -254.3037 -1445.3037
[29,] -4284.5037 -254.3037
[30,] -9687.3037 -4284.5037
[31,] -10046.3037 -9687.3037
[32,] -11495.5037 -10046.3037
[33,] -21621.4074 -11495.5037
[34,] -19588.0074 -21621.4074
[35,] -19322.6074 -19588.0074
[36,] -21349.7037 -19322.6074
[37,] -19897.9037 -21349.7037
[38,] -21558.9037 -19897.9037
[39,] -20201.3037 -21558.9037
[40,] -20747.3037 -20201.3037
[41,] -23421.5037 -20747.3037
[42,] -27805.3037 -23421.5037
[43,] -26314.3037 -27805.3037
[44,] -31998.5037 -26314.3037
[45,] -22924.8889 -31998.5037
[46,] -22980.4889 -22924.8889
[47,] -24753.0889 -22980.4889
[48,] -21814.1852 -24753.0889
[49,] -18213.3852 -21814.1852
[50,] -11488.3852 -18213.3852
[51,] -4311.7852 -11488.3852
[52,] 2260.2148 -4311.7852
[53,] 7116.0148 2260.2148
[54,] 15390.2148 7116.0148
[55,] 17793.2148 15390.2148
[56,] 17107.0148 17793.2148
[57,] 22229.1111 17107.0148
[58,] 21795.5111 22229.1111
[59,] 22794.9111 21795.5111
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 12501.0963 10589.2963
2 11706.0963 12501.0963
3 10608.6963 11706.0963
4 7839.6963 10608.6963
5 8642.4963 7839.6963
6 10151.6963 8642.4963
7 6750.6963 10151.6963
8 9587.4963 6750.6963
9 7999.5926 9587.4963
10 9114.9926 7999.5926
11 10870.3926 9114.9926
12 16287.2963 10870.3926
13 17645.0963 16287.2963
14 16496.0963 17645.0963
15 15349.6963 16496.0963
16 10901.6963 15349.6963
17 11947.4963 10901.6963
18 11950.6963 11947.4963
19 11816.6963 11950.6963
20 16799.4963 11816.6963
21 14317.5926 16799.4963
22 11657.9926 14317.5926
23 10410.3926 11657.9926
24 16287.2963 10410.3926
25 7965.0963 16287.2963
26 4845.0963 7965.0963
27 -1445.3037 4845.0963
28 -254.3037 -1445.3037
29 -4284.5037 -254.3037
30 -9687.3037 -4284.5037
31 -10046.3037 -9687.3037
32 -11495.5037 -10046.3037
33 -21621.4074 -11495.5037
34 -19588.0074 -21621.4074
35 -19322.6074 -19588.0074
36 -21349.7037 -19322.6074
37 -19897.9037 -21349.7037
38 -21558.9037 -19897.9037
39 -20201.3037 -21558.9037
40 -20747.3037 -20201.3037
41 -23421.5037 -20747.3037
42 -27805.3037 -23421.5037
43 -26314.3037 -27805.3037
44 -31998.5037 -26314.3037
45 -22924.8889 -31998.5037
46 -22980.4889 -22924.8889
47 -24753.0889 -22980.4889
48 -21814.1852 -24753.0889
49 -18213.3852 -21814.1852
50 -11488.3852 -18213.3852
51 -4311.7852 -11488.3852
52 2260.2148 -4311.7852
53 7116.0148 2260.2148
54 15390.2148 7116.0148
55 17793.2148 15390.2148
56 17107.0148 17793.2148
57 22229.1111 17107.0148
58 21795.5111 22229.1111
59 22794.9111 21795.5111
> 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/72vv71258742239.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/8pvzj1258742239.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/9n0ki1258742239.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/108bvo1258742239.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/11d9o91258742239.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/12zslu1258742239.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/1338nr1258742239.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/14uyb91258742239.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/156p0m1258742239.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/16jb4o1258742239.tab")
+ }
>
> system("convert tmp/1hje01258742239.ps tmp/1hje01258742239.png")
> system("convert tmp/2l78t1258742239.ps tmp/2l78t1258742239.png")
> system("convert tmp/3vo3s1258742239.ps tmp/3vo3s1258742239.png")
> system("convert tmp/4j86l1258742239.ps tmp/4j86l1258742239.png")
> system("convert tmp/5u8co1258742239.ps tmp/5u8co1258742239.png")
> system("convert tmp/6vwoe1258742239.ps tmp/6vwoe1258742239.png")
> system("convert tmp/72vv71258742239.ps tmp/72vv71258742239.png")
> system("convert tmp/8pvzj1258742239.ps tmp/8pvzj1258742239.png")
> system("convert tmp/9n0ki1258742239.ps tmp/9n0ki1258742239.png")
> system("convert tmp/108bvo1258742239.ps tmp/108bvo1258742239.png")
>
>
> proc.time()
user system elapsed
2.414 1.528 2.805