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(416.25,1111.92,398.35,1131.13,400.00,1144.94,427.25,1113.89,391.25,1107.30,397.20,1120.68,394.80,1140.84,391.50,1101.72,407.65,1104.24,418.10,1114.58,429.10,1130.20,452.85,1173.78,427.75,1211.92,420.90,1181.27,433.45,1203.60,427.15,1180.59,427.90,1156.85,415.35,1191.50,432.60,1191.33,431.65,1234.18,439.60,1220.33,466.10,1228.81,459.50,1207.01,499.75,1249.48,530.00,1248.29,568.25,1280.08,564.25,1280.66,587.00,1302.88,661.00,1310.61,625.00,1270.05,622.95,1270.06,637.25,1278.53,621.05,1303.80,600.60,1335.83,614.10,1377.76,648.75,1400.63,639.75,1418.03,660.20,1437.90,670.40,1406.80,658.25,1420.83,673.60,1482.37,666.50,1530.63,654.75,1504.66,665.75,1455.18,672.00,1473.96,742.50,1527.29,790.25,1545.79,784.25,1479.63,846.75,1467.97,914.75,1378.60,988.50,1330.45,887.75,1326.41,853.00,1385.97,888.25,1399.62,937.50,1276.69,912.50,1269.42,822.25,1287.83,880.00,1164.17,729.50,968.67,778.00,888.61),dim=c(2,60),dimnames=list(c('Gold','S&P500'),1:60))
> y <- array(NA,dim=c(2,60),dimnames=list(c('Gold','S&P500'),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 = '2'
> #'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
S&P500 Gold M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 1111.92 416.25 1 0 0 0 0 0 0 0 0 0 0 1
2 1131.13 398.35 0 1 0 0 0 0 0 0 0 0 0 2
3 1144.94 400.00 0 0 1 0 0 0 0 0 0 0 0 3
4 1113.89 427.25 0 0 0 1 0 0 0 0 0 0 0 4
5 1107.30 391.25 0 0 0 0 1 0 0 0 0 0 0 5
6 1120.68 397.20 0 0 0 0 0 1 0 0 0 0 0 6
7 1140.84 394.80 0 0 0 0 0 0 1 0 0 0 0 7
8 1101.72 391.50 0 0 0 0 0 0 0 1 0 0 0 8
9 1104.24 407.65 0 0 0 0 0 0 0 0 1 0 0 9
10 1114.58 418.10 0 0 0 0 0 0 0 0 0 1 0 10
11 1130.20 429.10 0 0 0 0 0 0 0 0 0 0 1 11
12 1173.78 452.85 0 0 0 0 0 0 0 0 0 0 0 12
13 1211.92 427.75 1 0 0 0 0 0 0 0 0 0 0 13
14 1181.27 420.90 0 1 0 0 0 0 0 0 0 0 0 14
15 1203.60 433.45 0 0 1 0 0 0 0 0 0 0 0 15
16 1180.59 427.15 0 0 0 1 0 0 0 0 0 0 0 16
17 1156.85 427.90 0 0 0 0 1 0 0 0 0 0 0 17
18 1191.50 415.35 0 0 0 0 0 1 0 0 0 0 0 18
19 1191.33 432.60 0 0 0 0 0 0 1 0 0 0 0 19
20 1234.18 431.65 0 0 0 0 0 0 0 1 0 0 0 20
21 1220.33 439.60 0 0 0 0 0 0 0 0 1 0 0 21
22 1228.81 466.10 0 0 0 0 0 0 0 0 0 1 0 22
23 1207.01 459.50 0 0 0 0 0 0 0 0 0 0 1 23
24 1249.48 499.75 0 0 0 0 0 0 0 0 0 0 0 24
25 1248.29 530.00 1 0 0 0 0 0 0 0 0 0 0 25
26 1280.08 568.25 0 1 0 0 0 0 0 0 0 0 0 26
27 1280.66 564.25 0 0 1 0 0 0 0 0 0 0 0 27
28 1302.88 587.00 0 0 0 1 0 0 0 0 0 0 0 28
29 1310.61 661.00 0 0 0 0 1 0 0 0 0 0 0 29
30 1270.05 625.00 0 0 0 0 0 1 0 0 0 0 0 30
31 1270.06 622.95 0 0 0 0 0 0 1 0 0 0 0 31
32 1278.53 637.25 0 0 0 0 0 0 0 1 0 0 0 32
33 1303.80 621.05 0 0 0 0 0 0 0 0 1 0 0 33
34 1335.83 600.60 0 0 0 0 0 0 0 0 0 1 0 34
35 1377.76 614.10 0 0 0 0 0 0 0 0 0 0 1 35
36 1400.63 648.75 0 0 0 0 0 0 0 0 0 0 0 36
37 1418.03 639.75 1 0 0 0 0 0 0 0 0 0 0 37
38 1437.90 660.20 0 1 0 0 0 0 0 0 0 0 0 38
39 1406.80 670.40 0 0 1 0 0 0 0 0 0 0 0 39
40 1420.83 658.25 0 0 0 1 0 0 0 0 0 0 0 40
41 1482.37 673.60 0 0 0 0 1 0 0 0 0 0 0 41
42 1530.63 666.50 0 0 0 0 0 1 0 0 0 0 0 42
43 1504.66 654.75 0 0 0 0 0 0 1 0 0 0 0 43
44 1455.18 665.75 0 0 0 0 0 0 0 1 0 0 0 44
45 1473.96 672.00 0 0 0 0 0 0 0 0 1 0 0 45
46 1527.29 742.50 0 0 0 0 0 0 0 0 0 1 0 46
47 1545.79 790.25 0 0 0 0 0 0 0 0 0 0 1 47
48 1479.63 784.25 0 0 0 0 0 0 0 0 0 0 0 48
49 1467.97 846.75 1 0 0 0 0 0 0 0 0 0 0 49
50 1378.60 914.75 0 1 0 0 0 0 0 0 0 0 0 50
51 1330.45 988.50 0 0 1 0 0 0 0 0 0 0 0 51
52 1326.41 887.75 0 0 0 1 0 0 0 0 0 0 0 52
53 1385.97 853.00 0 0 0 0 1 0 0 0 0 0 0 53
54 1399.62 888.25 0 0 0 0 0 1 0 0 0 0 0 54
55 1276.69 937.50 0 0 0 0 0 0 1 0 0 0 0 55
56 1269.42 912.50 0 0 0 0 0 0 0 1 0 0 0 56
57 1287.83 822.25 0 0 0 0 0 0 0 0 1 0 0 57
58 1164.17 880.00 0 0 0 0 0 0 0 0 0 1 0 58
59 968.67 729.50 0 0 0 0 0 0 0 0 0 0 1 59
60 888.61 778.00 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) Gold M1 M2 M3 M4
1.050e+03 8.864e-02 9.892e+01 8.361e+01 6.977e+01 6.296e+01
M5 M6 M7 M8 M9 M10
7.865e+01 8.911e+01 5.877e+01 4.627e+01 5.417e+01 4.404e+01
M11 t
1.363e+01 3.668e+00
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-450.72 -52.32 -11.01 64.70 239.43
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.050e+03 1.170e+02 8.981 1.11e-11 ***
Gold 8.864e-02 3.301e-01 0.269 0.789
M1 9.892e+01 8.745e+01 1.131 0.264
M2 8.361e+01 8.801e+01 0.950 0.347
M3 6.977e+01 8.859e+01 0.788 0.435
M4 6.296e+01 8.689e+01 0.725 0.472
M5 7.865e+01 8.651e+01 0.909 0.368
M6 8.911e+01 8.594e+01 1.037 0.305
M7 5.877e+01 8.588e+01 0.684 0.497
M8 4.627e+01 8.557e+01 0.541 0.591
M9 5.417e+01 8.548e+01 0.634 0.529
M10 4.404e+01 8.542e+01 0.516 0.609
M11 1.363e+01 8.557e+01 0.159 0.874
t 3.668e+00 3.375e+00 1.087 0.283
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 134.9 on 46 degrees of freedom
Multiple R-squared: 0.3085, Adjusted R-squared: 0.1131
F-statistic: 1.579 on 13 and 46 DF, p-value: 0.1266
> 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.090578e-03 2.181156e-03 0.9989094
[2,] 8.589234e-05 1.717847e-04 0.9999141
[3,] 5.976529e-06 1.195306e-05 0.9999940
[4,] 4.563929e-05 9.127858e-05 0.9999544
[5,] 1.377426e-05 2.754853e-05 0.9999862
[6,] 2.715360e-06 5.430720e-06 0.9999973
[7,] 3.797805e-07 7.595609e-07 0.9999996
[8,] 5.128406e-08 1.025681e-07 0.9999999
[9,] 1.799961e-08 3.599922e-08 1.0000000
[10,] 3.222841e-09 6.445683e-09 1.0000000
[11,] 4.337971e-10 8.675941e-10 1.0000000
[12,] 1.690469e-10 3.380938e-10 1.0000000
[13,] 6.981484e-11 1.396297e-10 1.0000000
[14,] 6.162753e-11 1.232551e-10 1.0000000
[15,] 4.950746e-11 9.901492e-11 1.0000000
[16,] 4.304762e-11 8.609523e-11 1.0000000
[17,] 1.763448e-10 3.526896e-10 1.0000000
[18,] 9.584509e-10 1.916902e-09 1.0000000
[19,] 1.118944e-07 2.237889e-07 0.9999999
[20,] 7.868815e-06 1.573763e-05 0.9999921
[21,] 3.536114e-04 7.072229e-04 0.9996464
[22,] 4.570797e-04 9.141594e-04 0.9995429
[23,] 1.947091e-04 3.894182e-04 0.9998053
[24,] 1.572335e-04 3.144670e-04 0.9998428
[25,] 8.592915e-03 1.718583e-02 0.9914071
[26,] 5.250988e-02 1.050198e-01 0.9474901
[27,] 1.122097e-01 2.244193e-01 0.8877903
> postscript(file="/var/www/html/rcomp/tmp/1mo341259331723.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/2z7wk1259331723.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/3qyvy1259331723.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/4f8f01259331723.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/5c4yp1259331723.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
-77.8689182 -45.4347593 -21.5958643 -51.9182249 -74.6740093 -75.9536177
7 8 9 10 11 12
-28.9091199 -58.8966248 -69.3833708 -53.5074664 -12.1168246 39.3202982
13 14 15 16 17 18
-22.8993158 -41.3046674 -9.9119863 -29.2203582 -72.3837905 -50.7534945
19 20 21 22 23 24
-25.7808411 25.9933419 -0.1365276 12.4566489 17.9874160 66.8519214
25 26 27 28 29 30
-39.6041094 0.4327180 11.5424487 34.8889837 16.7024163 -34.7985980
31 32 33 34 35 36
-7.9351254 8.1072444 23.2381150 63.5431030 131.0221361 160.7830450
37 38 39 40 41 42
76.3962709 106.0909523 84.2619455 102.5121380 143.3345109 178.0916999
43 44 45 46 47 48
179.8350143 138.2199076 144.8707321 198.4135954 239.4265892 183.7608556
49 50 51 52 53 54
63.9760725 -19.7842436 -64.2965437 -56.2625386 -12.9791274 -16.5859898
55 56 57 58 59 60
-117.2099278 -113.4238691 -98.5889487 -220.9058810 -376.3193168 -450.7161202
> postscript(file="/var/www/html/rcomp/tmp/6l3cq1259331723.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 -77.8689182 NA
1 -45.4347593 -77.8689182
2 -21.5958643 -45.4347593
3 -51.9182249 -21.5958643
4 -74.6740093 -51.9182249
5 -75.9536177 -74.6740093
6 -28.9091199 -75.9536177
7 -58.8966248 -28.9091199
8 -69.3833708 -58.8966248
9 -53.5074664 -69.3833708
10 -12.1168246 -53.5074664
11 39.3202982 -12.1168246
12 -22.8993158 39.3202982
13 -41.3046674 -22.8993158
14 -9.9119863 -41.3046674
15 -29.2203582 -9.9119863
16 -72.3837905 -29.2203582
17 -50.7534945 -72.3837905
18 -25.7808411 -50.7534945
19 25.9933419 -25.7808411
20 -0.1365276 25.9933419
21 12.4566489 -0.1365276
22 17.9874160 12.4566489
23 66.8519214 17.9874160
24 -39.6041094 66.8519214
25 0.4327180 -39.6041094
26 11.5424487 0.4327180
27 34.8889837 11.5424487
28 16.7024163 34.8889837
29 -34.7985980 16.7024163
30 -7.9351254 -34.7985980
31 8.1072444 -7.9351254
32 23.2381150 8.1072444
33 63.5431030 23.2381150
34 131.0221361 63.5431030
35 160.7830450 131.0221361
36 76.3962709 160.7830450
37 106.0909523 76.3962709
38 84.2619455 106.0909523
39 102.5121380 84.2619455
40 143.3345109 102.5121380
41 178.0916999 143.3345109
42 179.8350143 178.0916999
43 138.2199076 179.8350143
44 144.8707321 138.2199076
45 198.4135954 144.8707321
46 239.4265892 198.4135954
47 183.7608556 239.4265892
48 63.9760725 183.7608556
49 -19.7842436 63.9760725
50 -64.2965437 -19.7842436
51 -56.2625386 -64.2965437
52 -12.9791274 -56.2625386
53 -16.5859898 -12.9791274
54 -117.2099278 -16.5859898
55 -113.4238691 -117.2099278
56 -98.5889487 -113.4238691
57 -220.9058810 -98.5889487
58 -376.3193168 -220.9058810
59 -450.7161202 -376.3193168
60 NA -450.7161202
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -45.4347593 -77.8689182
[2,] -21.5958643 -45.4347593
[3,] -51.9182249 -21.5958643
[4,] -74.6740093 -51.9182249
[5,] -75.9536177 -74.6740093
[6,] -28.9091199 -75.9536177
[7,] -58.8966248 -28.9091199
[8,] -69.3833708 -58.8966248
[9,] -53.5074664 -69.3833708
[10,] -12.1168246 -53.5074664
[11,] 39.3202982 -12.1168246
[12,] -22.8993158 39.3202982
[13,] -41.3046674 -22.8993158
[14,] -9.9119863 -41.3046674
[15,] -29.2203582 -9.9119863
[16,] -72.3837905 -29.2203582
[17,] -50.7534945 -72.3837905
[18,] -25.7808411 -50.7534945
[19,] 25.9933419 -25.7808411
[20,] -0.1365276 25.9933419
[21,] 12.4566489 -0.1365276
[22,] 17.9874160 12.4566489
[23,] 66.8519214 17.9874160
[24,] -39.6041094 66.8519214
[25,] 0.4327180 -39.6041094
[26,] 11.5424487 0.4327180
[27,] 34.8889837 11.5424487
[28,] 16.7024163 34.8889837
[29,] -34.7985980 16.7024163
[30,] -7.9351254 -34.7985980
[31,] 8.1072444 -7.9351254
[32,] 23.2381150 8.1072444
[33,] 63.5431030 23.2381150
[34,] 131.0221361 63.5431030
[35,] 160.7830450 131.0221361
[36,] 76.3962709 160.7830450
[37,] 106.0909523 76.3962709
[38,] 84.2619455 106.0909523
[39,] 102.5121380 84.2619455
[40,] 143.3345109 102.5121380
[41,] 178.0916999 143.3345109
[42,] 179.8350143 178.0916999
[43,] 138.2199076 179.8350143
[44,] 144.8707321 138.2199076
[45,] 198.4135954 144.8707321
[46,] 239.4265892 198.4135954
[47,] 183.7608556 239.4265892
[48,] 63.9760725 183.7608556
[49,] -19.7842436 63.9760725
[50,] -64.2965437 -19.7842436
[51,] -56.2625386 -64.2965437
[52,] -12.9791274 -56.2625386
[53,] -16.5859898 -12.9791274
[54,] -117.2099278 -16.5859898
[55,] -113.4238691 -117.2099278
[56,] -98.5889487 -113.4238691
[57,] -220.9058810 -98.5889487
[58,] -376.3193168 -220.9058810
[59,] -450.7161202 -376.3193168
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -45.4347593 -77.8689182
2 -21.5958643 -45.4347593
3 -51.9182249 -21.5958643
4 -74.6740093 -51.9182249
5 -75.9536177 -74.6740093
6 -28.9091199 -75.9536177
7 -58.8966248 -28.9091199
8 -69.3833708 -58.8966248
9 -53.5074664 -69.3833708
10 -12.1168246 -53.5074664
11 39.3202982 -12.1168246
12 -22.8993158 39.3202982
13 -41.3046674 -22.8993158
14 -9.9119863 -41.3046674
15 -29.2203582 -9.9119863
16 -72.3837905 -29.2203582
17 -50.7534945 -72.3837905
18 -25.7808411 -50.7534945
19 25.9933419 -25.7808411
20 -0.1365276 25.9933419
21 12.4566489 -0.1365276
22 17.9874160 12.4566489
23 66.8519214 17.9874160
24 -39.6041094 66.8519214
25 0.4327180 -39.6041094
26 11.5424487 0.4327180
27 34.8889837 11.5424487
28 16.7024163 34.8889837
29 -34.7985980 16.7024163
30 -7.9351254 -34.7985980
31 8.1072444 -7.9351254
32 23.2381150 8.1072444
33 63.5431030 23.2381150
34 131.0221361 63.5431030
35 160.7830450 131.0221361
36 76.3962709 160.7830450
37 106.0909523 76.3962709
38 84.2619455 106.0909523
39 102.5121380 84.2619455
40 143.3345109 102.5121380
41 178.0916999 143.3345109
42 179.8350143 178.0916999
43 138.2199076 179.8350143
44 144.8707321 138.2199076
45 198.4135954 144.8707321
46 239.4265892 198.4135954
47 183.7608556 239.4265892
48 63.9760725 183.7608556
49 -19.7842436 63.9760725
50 -64.2965437 -19.7842436
51 -56.2625386 -64.2965437
52 -12.9791274 -56.2625386
53 -16.5859898 -12.9791274
54 -117.2099278 -16.5859898
55 -113.4238691 -117.2099278
56 -98.5889487 -113.4238691
57 -220.9058810 -98.5889487
58 -376.3193168 -220.9058810
59 -450.7161202 -376.3193168
> 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/7p9cc1259331723.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/89mxz1259331723.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/96ht41259331723.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/10fa3k1259331723.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/11ksix1259331723.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/12d2ba1259331723.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/139m071259331723.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/14mks41259331723.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/150mjd1259331723.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/164z8g1259331723.tab")
+ }
>
> system("convert tmp/1mo341259331723.ps tmp/1mo341259331723.png")
> system("convert tmp/2z7wk1259331723.ps tmp/2z7wk1259331723.png")
> system("convert tmp/3qyvy1259331723.ps tmp/3qyvy1259331723.png")
> system("convert tmp/4f8f01259331723.ps tmp/4f8f01259331723.png")
> system("convert tmp/5c4yp1259331723.ps tmp/5c4yp1259331723.png")
> system("convert tmp/6l3cq1259331723.ps tmp/6l3cq1259331723.png")
> system("convert tmp/7p9cc1259331723.ps tmp/7p9cc1259331723.png")
> system("convert tmp/89mxz1259331723.ps tmp/89mxz1259331723.png")
> system("convert tmp/96ht41259331723.ps tmp/96ht41259331723.png")
> system("convert tmp/10fa3k1259331723.ps tmp/10fa3k1259331723.png")
>
>
> proc.time()
user system elapsed
2.388 1.569 2.838