** R output *   Output edited.
*************
**  Problem A **

> setwd("C:/teaching/bs41202")
> da=read.table("m-jnj6105.txt")
> jnj=log(da[,2]+1)
> library(fSeries)
> basicStats(jnj)
                    Value
nobs        540.000000000
NAs           0.000000000
Minimum      -0.189442853
Maximum       0.192453389
1. Quartile  -0.028394337
3. Quartile   0.055451205
Mean          0.012186223
Median        0.013108700
Sum           6.580560235
SE Mean       0.002751256
LCL Mean      0.006781725
UCL Mean      0.017590720
Variance      0.004087480
Stdev         0.063933403
Skewness     -0.119286294
Kurtosis      0.084330569
> Box.test(jnj,lag=10,type='Ljung')
        Box-Ljung test
data:  jnj 
X-squared = 10.2624, df = 10, p-value = 0.4178

** Problem B **
> da=read.table("gas-sent.txt")
> y=da[,1]
> x=as.matrix(da[,2:5])

> m1=lm(y~x)
> summary(m1)
lm(formula = y ~ x)

Residuals:
    Min      1Q  Median      3Q     Max 
-36.521 -13.817  -2.781  13.556  35.832 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   1.8581    11.5561   0.161  0.87282    
xV2           5.8742     0.2905  20.219  < 2e-16 ***
xV3           1.4052     0.2928   4.799 1.16e-05 ***
xV4           1.3154     0.5787   2.273  0.02675 *  
xV5         -15.8571     5.3344  -2.973  0.00429 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Residual standard error: 18.32 on 58 degrees of freedom
Multiple R-Squared: 0.9521,     Adjusted R-squared: 0.9488 
F-statistic: 288.1 on 4 and 58 DF,  p-value: < 2.2e-16 

> a1=acf(m1$residuals,lag=10)
> a1$acf
 [1,]  1.00000000
 [2,]  0.51514249
 [3,]  0.27614439
 [4,]  0.25870432
 [5,]  0.27702979
 [6,]  0.10484658
 [7,]  0.26511473
 [8,]  0.35066692
  ..
[11,] -0.01653887

> a2=pacf(m1$residuals,lag.max=10)
> a2$acf
 [1,]  0.51514249
 [2,]  0.01466401
 [3,]  0.15110560
 [4,]  0.11923855
 [5,] -0.14493573
 [6,]  0.32408261
 [7,]  0.10540219
 [8,] -0.29272433
 [9,] -0.02825233
[10,] -0.10040808

> m2=arima(y,xreg=x,order=c(1,0,0),seasonal=list(order=c(1,0,0),period=7))
> m2
arima(x=y,order=c(1,0,0),seasonal=list(order=c(1,0,0),period = 7), 
    xreg = x)

Coefficients:
         ar1    sar1  intercept      V2      V3      V4        V5
      0.5414  0.3753     7.7211  5.7066  1.4087  1.2008  -10.0592
s.e.  0.1067  0.1157    12.5670  0.2195  0.2211  0.3840    6.0951

sigma^2 estimated as 188.2:  log likelihood = -255.09,  aic = 526.18

> Box.test(m2$residuals,lag=14,type='Ljung')
        Box-Ljung test
data:  m2$residuals 
X-squared = 8.1352, df = 14, p-value = 0.8822

** Problem C **
> source("garchoxfit.R")
> da=read.table("d-qcomsp6205.txt")
> qcom=log(da[,2]+1)
> m1=garchOxFit(formula.mean=~arma(0,0),formula.var=~garch(1,1),series=qcom,cond.dist="t")
 ********************
 ** SPECIFICATIONS **
 ********************
Mean Equation : ARMA (0, 0) model.
No regressor in the mean
Variance Equation : GARCH (1, 1) model.
 No regressor in the variance
The distribution is a Student distribution, with 6.00115 degrees of freedom.

Strong convergence using numerical derivatives
Log-likelihood = 6894.7

 Maximum Likelihood Estimation (Std.Errors based on Numerical OPG matrix)
                  Coefficient  Std.Error  t-value  t-prob
Cst(M)               0.000337 0.00048545   0.6936  0.4880
Cst(V)               0.000000  0.0062740     0.00  1.0000
ARCH(Alpha1)         0.017537  0.0028534    6.146  0.0000
GARCH(Beta1)         0.982408  0.0027840    352.9  0.0000
Student(DF)          6.001147    0.52359    11.46  0.0000

No. Observations :      3530  No. Parameters  :         5
Mean (Y)         :   0.00116  Variance (Y)    :   0.00150
Skewness (Y)     :   0.49689  Kurtosis (Y)    :   7.35492
Log Likelihood   :  6894.700  Alpha[1]+Beta[1]:   0.99995

Warning : To avoid numerical problems, the estimated parameter
Cst(V), and its std.Error have been multiplied by 10^4.

*** GJR Model
> m2=garchOxFit(formula.mean=~arma(0,0),formula.var=~gjr(1,1),series=qcom,cond.dist="t")
 ********************
 ** SPECIFICATIONS **
 ********************
Mean Equation : ARMA (0, 0) model.
No regressor in the mean
Variance Equation : GJR (1, 1) model.
 No regressor in the variance
The distribution is a Student distribution, with 6.28931 degrees of freedom.

Strong convergence using numerical derivatives
Log-likelihood = 6899.56

 Maximum Likelihood Estimation (Std.Errors based on Numerical OPG matrix)
                  Coefficient  Std.Error  t-value  t-prob
Cst(M)               0.000318 0.00048983   0.6502  0.5156
Cst(V)               0.048013   0.021026    2.283  0.0225
ARCH(Alpha1)         0.031232  0.0061617    5.069  0.0000
GARCH(Beta1)         0.947233  0.0064067    147.8  0.0000
GJR(Gamma1)          0.044792   0.011587    3.866  0.0001
Student(DF)          6.289314    0.55337    11.37  0.0000

No. Observations :      3530  No. Parameters  :         6
Mean (Y)         :   0.00116  Variance (Y)    :   0.00150
Skewness (Y)     :   0.49689  Kurtosis (Y)    :   7.35492
Log Likelihood   :  6899.564

Warning : To avoid numerical problems, the estimated parameter
Cst(V), and its std.Error have been multiplied by 10^4.

***************
** FORECASTS **
***************
Number of Forecasts: 15

Horizon       Mean   Variance
      1  0.0003185  0.0003102
      2  0.0003185  0.0003105
      3  0.0003185  0.0003107
      4  0.0003185   0.000311
   .....
     15  0.0003185   0.000314
***********
** TESTS **
***********
                   Statistic       t-Test      P-Value
Skewness             0.51209       12.426  1.8812e-035
Excess Kurtosis       3.4910       42.368      0.00000
Jarque-Bera           1946.8         .NaN      0.00000
---------------
Information Criterium (to be minimized)
Akaike         -3.905702  Shibata        -3.905708
Schwarz         -3.895217  Hannan-Quinn    -3.901962
---------------
Q-Statistics on Standardized Residuals
  Q( 10) =  9.64676   [0.4720109]
  Q( 15) =  13.9482   [0.5294635]
  Q( 20) =  22.5577   [0.3110326]
--------------
Q-Statistics on Squared Standardized Residuals
  --> P-values adjusted by 2 degree(s) of freedom 
  Q( 10) =  12.0071   [0.1508875]
  Q( 15) =  15.2297   [0.2932329]
  Q( 20) =  22.4805   [0.2113499]

*** EGARCH
> m3=garchOxFit(formula.mean=~arma(0,0),formula.var=~egarch(1,1),series=qcom,cond.dist="t")
 ********************
 ** SPECIFICATIONS **
 ********************
Mean Equation : ARMA (0, 0) model.
No regressor in the mean
Variance Equation : EGARCH (1, 1) model.
 No regressor in the variance
The distribution is a Student distribution, with 4.37418 degrees of freedom.

Strong convergence using numerical derivatives
Log-likelihood = 6855.83

 Maximum Likelihood Estimation (Std.Errors based on Numerical OPG matrix)
                  Coefficient  Std.Error  t-value  t-prob
Cst(M)              -0.000216 0.00050934  -0.4234  0.6720
Cst(V)             100.000000     10686. 0.009358  0.9925
ARCH(Alpha1)         0.063729    0.19961   0.3193  0.7495
GARCH(Beta1)         0.884271   0.020319    43.52  0.0000
EGARCH(Theta1)      -0.092714   0.023598   -3.929  0.0001
EGARCH(Theta2)       0.262061   0.049642    5.279  0.0000
Student(DF)          4.374182    0.38194    11.45  0.0000

Warning : To avoid numerical problems, the estimated parameter
Cst(V), and its std.Error have been multiplied by 10^4.

***********
** TESTS **
***********
                   Statistic       t-Test      P-Value
Skewness             0.46134       11.195  4.3270e-029
Excess Kurtosis       3.2432       39.361      0.00000
Jarque-Bera           1672.3         .NaN      0.00000
---------------
Information Criterium (to be minimized)
Akaike         -3.880357  Shibata        -3.880364
Schwarz         -3.868123  Hannan-Quinn    -3.875993
---------------
Q-Statistics on Standardized Residuals
  Q( 10) =  9.22525   [0.5108680]
  Q( 15) =  15.2923   [0.4305774]
  Q( 20) =  24.5423   [0.2195069]
--------------
Q-Statistics on Squared Standardized Residuals
  --> P-values adjusted by 2 degree(s) of freedom 
  Q( 10) =  7.15925   [0.5195528]
  Q( 15) =  13.7496   [0.3917142]
  Q( 20) =  42.1940   [0.0010389]
--------------

** Problem D **
> da=read.table("q-earn-cat.txt")
> cat=da[,4]
> plot(cat,type='l')
> a1=acf(cat,lag.max=12)
> a1$acf
 [1,]  1.00000000
 [2,]  0.75849100
 [3,]  0.65390191
 [4,]  0.50004431
 [5,]  0.46466130
 [6,]  0.29776771
 [7,]  0.23462919
 [8,]  0.09848228
   ...
[13,] -0.08300225

> a2=acf(diff(cat),lag.max=12)
> a2$acf
 [1,]  1.0000000
 [2,] -0.5487003
 [3,]  0.4710748
 [4,] -0.4357344
 [5,]  0.5822293
 [6,] -0.4716251
 [7,]  0.4714451
 [8,] -0.4427083
 [9,]  0.4938262
[10,] -0.4681817
[11,]  0.4049203
[12,] -0.3114033
[13,]  0.3160038

> a3=acf(diff(diff(cat),4),lag.max=12)
> a3$acf
 [1,]  1.000000000
 [2,] -0.270698998
 [3,] -0.014342673
 [4,]  0.200337950
 [5,] -0.379147506
 [6,]  0.106760619
 [7,]  0.064983202
 [8,] -0.099839183
 [9,]  0.010919858
[10,] -0.082111896
[11,]  0.095345210
[12,] -0.002809605
[13,] -0.014018064

> m1=arima(cat,order=c(0,1,1),seasonal=list(order=c(0,1,1),period=4))
> m1
arima(x=cat,order=c(0,1,1),seasonal=list(order=c(0,1,1),period=4))

Coefficients:
          ma1     sma1
      -0.2158  -0.4736
s.e.   0.1319   0.1505

sigma^2 estimated as 0.008644:  log likelihood = 48.25,  aic = -90.5

> tsdiag(m1)  % Plots show the model fits the data well.

> Box.test(m1$residuals,lag=10,type='Ljung')
        Box-Ljung test
data:  m1$residuals 
X-squared = 4.1507, df = 10, p-value = 0.9403

> predict(m1,5)
$pred
Time Series:
Start = 57 
[1] 1.047657 1.043492 1.291565 1.163953 1.269883
$se
Time Series:
Start = 57 
[1] 0.0929724 0.1181501 0.1388345 0.1568139 0.1985889

*** S-Plus output
** S-Plus output ** Output is edited.

** Problem A **
> module(finmetrics)
> da=read.table("m-jnj6105.txt")

> jnj=log(da[,2]+1)  % log returns
> summaryStats(jnj)

Sample Quantiles:
     min       1Q  median      3Q    max
 -0.1894 -0.02839 0.01311 0.05545 0.1925

Sample Moments:
    mean     std skewness kurtosis
 0.01219 0.06393  -0.1196    3.096

> autocorTest(jnj,lag=10)
Test for Autocorrelation: Ljung-Box
 
Null Hypothesis: no autocorrelation
 
Test Statistics:
Test Stat 10.2624,  p.value  0.4178
 
** Problem B **
> da=read.table("gas.dat")
> y=da[,1]
> x=da[,2:5]
> x=as.matrix(x)
> m1=lm(y~x)
> summary(m1)
Call: lm(formula = y ~ x)
 
Coefficients:
               Value Std. Error  t value Pr(>|t|)
(Intercept)   1.8581  11.5561     0.1608   0.8728
        xV2   5.8742   0.2905    20.2194   0.0000
        xV3   1.4052   0.2928     4.7994   0.0000
        xV4   1.3154   0.5787     2.2730   0.0267
        xV5 -15.8571   5.3344    -2.9726   0.0043
 
Residual standard error: 18.32 on 58 degrees of freedom
Multiple R-Squared: 0.9521
F-statistic: 288.1 on 4 and 58 degrees of freedom, the p-value is 0

> autocorTest(m1$residuals,lag=10)
Test for Autocorrelation: Ljung-Box
 
Null Hypothesis: no autocorrelation
 
Test Statistics:
Test Stat 48.1743,  p.value  0.0000
 
Dist. under Null: chi-square with 10 degrees of freedom
   Total Observ.: 63

> acf(m1$residuals,lag=10)
Call: acf(x = m1$residuals, lag.max = 10)
Autocorrelation matrix:
   lag      m1
 1   0  1.0000
 2   1  0.5151
 3   2  0.2761
 4   3  0.2587
 5   4  0.2770
 6   5  0.1048
 7   6  0.2651
 8   7  0.3507
   ...
11  10 -0.0165
> acf(m1$residuals,lag=10,type='partial')
Call: acf(x = m1$residuals, lag.max = 10, type = "partial")
Partial Correlation matrix:
   lag      m1
 1   1  0.5151
 2   2  0.0147
 3   3  0.1511
 4   4  0.1192
 5   5 -0.1449
 6   6  0.3241
 7   7  0.1054
 8   8 -0.2927
 9   9 -0.0283
10  10 -0.1004
> x1=cbind(rep(1,63),x)
> m2=arima.mle(y,xreg=x1,model=list(list(order=c(1,0,0)),list(order=c(1,0,0),period=7)))
> m2
Call: arima.mle(x = y, model = list(list(order = c(1, 0, 0)), 
list(order = c(1,0,0),period = 7)), xreg = x1)
Method:  Maximum Likelihood
Model #1:
        Model : 1 0 0
        Coefficients:
                AR : 0.55224
Model #2:
        Model : 1 0 0
        Period :  7
        Coefficients:
                AR : 0.3676
 
Variance-Covariance Matrix:
              ar(1)         ar(7)
ar(1)  0.0126465419 -0.0003892573
ar(7) -0.0003892573  0.0157369048
Coeffficients for regressor(s): x1
[1]   8.26462   5.80983   1.31381   1.30916 -11.45429
 
Optimizer has  converged
AIC: 461.51913

> m2$sigma2
[1] 200.1034

** Problem C **
> da=read.table("d-qcomsp9205.txt")
> qcom=log(da[,2]+1)
> autocorTest(qcom,lag=10)
Test for Autocorrelation: Ljung-Box
 
Null Hypothesis: no autocorrelation
 
Test Statistics:
Test Stat 16.8375,  p.value  0.0780
 
> archTest(qcom,lag=10)
Test for ARCH Effects: LM Test
 
Null Hypothesis: no ARCH effects
 
Test Statistics:
Test Stat 154.1304,  p.value   0.0000
 
> m1=garch(qcom~1,~garch(1,1),cond.dist="t")
> summary(m1)
garch(formula.mean = qcom ~ 1, formula.var =  ~ garch(1, 1), cond.dist = "t")
 
Mean Equation: qcom ~ 1
 
Conditional Variance Equation:  ~ garch(1, 1)
 
Conditional Distribution:  t
 with estimated parameter 6.05477 and standard error 0.5334613
--------------------------------------------------------------
Estimated Coefficients:
--------------------------------------------------------------
             Value Std.Error t value  Pr(>|t|)
       C 5.568e-04 4.854e-04   1.147 2.514e-01
       A 1.917e-06 1.367e-06   1.402 1.611e-01
 ARCH(1) 3.558e-02 5.139e-03   6.923 5.249e-12
GARCH(1) 9.638e-01 4.939e-03 195.141 0.000e+00
--------------------------------------------------------------
AIC(5) = -13775.4, BIC(5) = -13744.56
 
Ljung-Box test for standardized residuals:
--------------------------------------------------------------
 Statistic P-value Chi^2-d.f.
     11.59  0.4791         12
 
Ljung-Box test for squared standardized residuals:
--------------------------------------------------------------
 Statistic P-value Chi^2-d.f.
     13.84  0.3108         12
 
> m2=garch(qcom~1,~tgarch(1,1),cond.dist="t",leverage=T)
Warning messages:
  The estimated asymptotic variance is not well-defined. in: garch(qcom ~ 1,
        ~ tgarch(1, 1), cond.dist = "t", leverage = T)
> summary(m2)
garch(formula.mean=qcom ~ 1,formula.var=~tgarch(1, 1), leverage = T,
        cond.dist = "t")
 
Mean Equation: qcom ~ 1
 
Conditional Variance Equation:  ~ tgarch(1, 1)
 
Conditional Distribution:  t
 with estimated parameter 6.311176 and standard error 0.555231
--------------------------------------------------------------
Estimated Coefficients:
--------------------------------------------------------------
             Value Std.Error t value  Pr(>|t|)
       C 3.145e-04 4.915e-04   0.640 5.222e-01
       A 7.995e-06 2.871e-06   2.784 5.394e-03
 ARCH(1) 3.897e-02 7.469e-03   5.217 1.923e-07
GARCH(1) 9.315e-01 7.855e-03 118.591 0.000e+00
GAMMA(1) 6.069e-02 1.416e-02   4.287 1.861e-05
--------------------------------------------------------------
AIC(6) = -13786.23, BIC(6) = -13749.22
 
Ljung-Box test for standardized residuals:
--------------------------------------------------------------
 Statistic P-value Chi^2-d.f.
     10.18  0.6006         12
 
Ljung-Box test for squared standardized residuals:
--------------------------------------------------------------
 Statistic P-value Chi^2-d.f.
     10.44  0.5775         12
 
> predict(m2,5)
$series.pred:
[1] 0.0003145462 0.0003145462 0.0003145462 0.0003145462 0.0003145462
 
$sigma.pred:
[1] 0.01788012 0.01810924 0.01833568 0.01855952 0.01878087
 
> m3=garch(qcom~1,~egarch(1,1),leverage=T,cond.dist="t")
> summary(m3)
garch(formula.mean =qcom~1,formula.var=~egarch(1,1),leverage = T,
        cond.dist = "t")
 
Mean Equation: qcom ~ 1
 
Conditional Variance Equation:  ~ egarch(1, 1)
 
Conditional Distribution:  t
 with estimated parameter 6.467035 and standard error 0.5832481
--------------------------------------------------------------
Estimated Coefficients:
--------------------------------------------------------------
              Value Std.Error   t value  Pr(>|t|)
       C  0.0000605 0.0004995    0.1211 9.036e-01
       A -0.2869671 0.0431429   -6.6515 3.349e-11
 ARCH(1)  0.1833269 0.0197742    9.2710 0.000e+00
GARCH(1)  0.9777422 0.0051546  189.6828 0.000e+00
  LEV(1) -0.2876733 0.0637396   -4.5133 6.592e-06
--------------------------------------------------------------
AIC(6) = -13808.9, BIC(6) = -13771.89
 
Ljung-Box test for standardized residuals:
--------------------------------------------------------------
 Statistic P-value Chi^2-d.f.
     9.676  0.6444         12
 
Ljung-Box test for squared standardized residuals:
--------------------------------------------------------------
 Statistic P-value Chi^2-d.f.
     10.23   0.596         12

** Problem D ** 
> da=read.table("q-earn-cat.txt")
> cat=da[,4]
> acf(cat,lag.max=12)
Call: acf(x = cat, lag.max = 12)
Autocorrelation matrix:
   lag     cat
 1   0  1.0000
 2   1  0.7585
 3   2  0.6539
 4   3  0.5000
 5   4  0.4647
 6   5  0.2978
 7   6  0.2346
 8   7  0.0985
 9   8  0.0794
  ...
13  12 -0.0830
> acf(diff(cat),lag.max=12)
Call: acf(x = diff(cat), lag.max = 12)
Autocorrelation matrix:
   lag     cat
 1   0  1.0000
 2   1 -0.5487
 3   2  0.4711
 4   3 -0.4357
 5   4  0.5822
 6   5 -0.4716
 7   6  0.4714
 8   7 -0.4427
 9   8  0.4938
10   9 -0.4682
11  10  0.4049
12  11 -0.3114
13  12  0.3160
> acf(diff(diff(cat),4),lag.max=12)
Call: acf(x = diff(diff(cat), 4), lag.max = 12)
Autocorrelation matrix:
   lag     cat
 1   0  1.0000
 2   1 -0.2707
 3   2 -0.0143
 4   3  0.2003
 5   4 -0.3791
 6   5  0.1068
 7   6  0.0650
 8   7 -0.0998
 9   8  0.0109
   ....
13  12 -0.0140
> m1=arima.mle(cat,model=list(list(order=c(0,1,1)),list(order=c(0,1,1),period=4)))
> summary(m1)
Call: arima.mle(x=cat,model=list(list(order=c(0,1,1)),list(order=c(0,1,1),period=4)))
Method:  Maximum Likelihood  with likelihood conditional on  5 observations
 
Multiplicative ARIMA model --
Model component  1
ARIMA order:  0 1 1
 
Model component  2
ARIMA order:  0 1 1
Period:  4
 
       Value Std. Error t-value
ma(1) 0.2158     0.1367   1.578
ma(4) 0.4736     0.1233   3.840
 
Variance-Covariance Matrix:
              ma(1)         ma(4)
ma(1)  0.0186963618 -0.0001623637
ma(4) -0.0001623637  0.0152111779
 
Estimated innovations variance:  0.0086
 
Optimizer has  converged
AIC: -92.4987

> arima.diag(m1)  % Plots show the model fits the data well.

> For your information, the Ljung-Box statistics for the 
> residuals given Q(10) = 4.15 with p-value 0.94 (from R output). 

> arima.forecast(cat,model=m1$model,5)
$mean:
[1] 1.047657 1.043492 1.291565 1.163953 1.269883
 
$std.err:
[1] 0.09297246 0.11815018 0.13883463 0.15681403 0.19858915