Load Required Libraries

# Load Library

library(fGarch)
## NOTE: Packages 'fBasics', 'timeDate', and 'timeSeries' are no longer
## attached to the search() path when 'fGarch' is attached.
## 
## If needed attach them yourself in your R script by e.g.,
##         require("timeSeries")
library(rugarch)
## Loading required package: parallel
## 
## Attaching package: 'rugarch'
## The following object is masked from 'package:stats':
## 
##     sigma
library(fitdistrplus)
## Loading required package: MASS
## Loading required package: survival
## 
## Attaching package: 'fitdistrplus'
## The following object is masked from 'package:rugarch':
## 
##     fitdist

Load Data

data = read.csv("DATA25.csv")

Select The Stocks Column

stock_data = data$STOCK; stock_data
##   [1] 3.74 3.81 3.98 4.21 4.50 4.81 5.14 5.50 5.77 5.97 4.64 6.49 6.73 6.97 7.26
##  [16] 7.50 7.67 7.87 8.10 8.32 8.52 8.68 8.77 8.78 8.71 8.59 8.30 7.95 7.66 7.38
##  [31] 7.48 6.81 6.61 6.45 6.29 6.16 6.05 6.29 6.23 6.08 5.62 6.10 6.07 5.35 5.71
##  [46] 5.53 5.35 5.20 4.66 5.17 5.16 5.74 5.41 5.53 5.67 5.79 5.87 6.01 6.16 6.18
##  [61] 6.03 5.84 5.72 5.29 5.20 5.19 5.19 5.24 5.40 5.32 5.16 5.04 4.91 4.91 4.67
##  [76] 4.65 4.68 4.69 4.59 4.53 4.53 4.63 4.95 5.20 5.61 6.24 6.89 7.40 7.79 7.98
##  [91] 8.15 8.33 8.40 8.36 8.21 8.13 7.84 7.20 6.76 6.43 6.26 6.30 6.43 6.48 6.50
## [106] 6.47 6.44 6.46 6.59 6.72 6.88 6.87 6.77 6.58 6.42 6.31 6.29 6.34 6.54

Model Fit

M1 = ugarchspec(variance.model = list(model = "eGARCH", garchOrder =c(1,2)),mean.model = list(armaOrder = c(1,2), include.mean = FALSE),distribution.model = "norm");M1
## 
## *---------------------------------*
## *       GARCH Model Spec          *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## ------------------------------------
## GARCH Model      : eGARCH(1,2)
## Variance Targeting   : FALSE 
## 
## Conditional Mean Dynamics
## ------------------------------------
## Mean Model       : ARFIMA(1,0,2)
## Include Mean     : FALSE 
## GARCH-in-Mean        : FALSE 
## 
## Conditional Distribution
## ------------------------------------
## Distribution :  norm 
## Includes Skew    :  FALSE 
## Includes Shape   :  FALSE 
## Includes Lambda  :  FALSE

Log Returns

log_returns = diff(log(stock_data)); log_returns = na.omit(log_returns); log_returns
##   [1]  0.018543578  0.043652630  0.056180828  0.066614749  0.066619687
##   [6]  0.066355995  0.067695013  0.047923988  0.034074847 -0.252032561
##  [11]  0.335548164  0.036312613  0.035040081  0.040764604  0.032523192
##  [16]  0.022413595  0.025741447  0.028805999  0.026798193  0.023754086
##  [21]  0.018605188  0.010315278  0.001139601 -0.008004617 -0.013873055
##  [26] -0.034343221 -0.043083586 -0.037159945 -0.037238345  0.013459153
##  [31] -0.093840672 -0.029808466 -0.024503523 -0.025119060 -0.020884293
##  [36] -0.018018506  0.038902799 -0.009584738 -0.024371637 -0.078673032
##  [41]  0.081957107 -0.004930166 -0.126262044  0.065122463 -0.032031208
##  [46] -0.033091255 -0.028437935 -0.109643177  0.103857240 -0.001936109
##  [51]  0.106522631 -0.059210117  0.021938723  0.025001302  0.020943174
##  [56]  0.013722342  0.023570115  0.024652029  0.003241494 -0.024571261
##  [61] -0.032016214 -0.020761991 -0.078150560 -0.017159620 -0.001924928
##  [66]  0.000000000  0.009587801  0.030077455 -0.014925650 -0.030536724
##  [71] -0.023530497 -0.026132140  0.000000000 -0.050114870 -0.004291852
##  [76]  0.006430890  0.002134473 -0.021552558 -0.013158085  0.000000000
##  [81]  0.021834929  0.066830708  0.049271049  0.075892094  0.106429463
##  [86]  0.099090903  0.071408915  0.051360860  0.024097552  0.021079516
##  [91]  0.021845529  0.008368250 -0.004773279 -0.018105504 -0.009792000
##  [96] -0.036322089 -0.085157808 -0.063058136 -0.050048352 -0.026794353
## [101]  0.006369448  0.020424905  0.007745972  0.003081667 -0.004626068
## [106] -0.004647568  0.003100778  0.019924031  0.019534806  0.023530497
## [111] -0.001454546 -0.014663019 -0.028466342 -0.024616628 -0.017282441
## [116] -0.003174606  0.007917698  0.031058397

Plot Log Returns

plot(log_returns, type = "l", main = "Log Returns of Stock Data", xlab = "Time", ylab = "Log Returns")

Fit the Model to Log Returns

fit = ugarchfit(spec = M1, data = log_returns); fit
## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : eGARCH(1,2)
## Mean Model   : ARFIMA(1,0,2)
## Distribution : norm 
## 
## Optimal Parameters
## ------------------------------------
##         Estimate  Std. Error  t value Pr(>|t|)
## ar1     0.252329    0.061107   4.1293 0.000036
## ma1     0.065529    0.029304   2.2362 0.025339
## ma2     0.071288    0.011035   6.4604 0.000000
## omega  -0.981820    0.529695  -1.8536 0.063803
## alpha1  0.308509    0.160107   1.9269 0.053993
## beta1   0.453908    0.112428   4.0373 0.000054
## beta2   0.368841    0.114145   3.2313 0.001232
## gamma1  1.193263    0.264535   4.5108 0.000006
## 
## Robust Standard Errors:
##         Estimate  Std. Error  t value Pr(>|t|)
## ar1     0.252329    0.111610   2.2608 0.023771
## ma1     0.065529    0.037495   1.7477 0.080521
## ma2     0.071288    0.006183  11.5292 0.000000
## omega  -0.981820    0.924261  -1.0623 0.288111
## alpha1  0.308509    0.218095   1.4146 0.157196
## beta1   0.453908    0.094155   4.8208 0.000001
## beta2   0.368841    0.121213   3.0429 0.002343
## gamma1  1.193263    0.306985   3.8870 0.000101
## 
## LogLikelihood : 208.7915 
## 
## Information Criteria
## ------------------------------------
##                     
## Akaike       -3.4032
## Bayes        -3.2154
## Shibata      -3.4117
## Hannan-Quinn -3.3270
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                          statistic p-value
## Lag[1]                       0.136  0.7123
## Lag[2*(p+q)+(p+q)-1][8]      4.108  0.7304
## Lag[4*(p+q)+(p+q)-1][14]     7.665  0.4319
## d.o.f=3
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                          statistic p-value
## Lag[1]                    0.006453  0.9360
## Lag[2*(p+q)+(p+q)-1][8]   0.538370  0.9954
## Lag[4*(p+q)+(p+q)-1][14]  0.984129  0.9997
## d.o.f=3
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[4]    0.1533 0.500 2.000  0.6954
## ARCH Lag[6]    0.2979 1.461 1.711  0.9457
## ARCH Lag[8]    0.5146 2.368 1.583  0.9807
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  2.1183
## Individual Statistics:             
## ar1    0.3294
## ma1    0.4962
## ma2    0.1858
## omega  0.2480
## alpha1 0.1131
## beta1  0.2948
## beta2  0.2981
## gamma1 0.1879
## 
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic:          1.89 2.11 2.59
## Individual Statistic:     0.35 0.47 0.75
## 
## Sign Bias Test
## ------------------------------------
##                    t-value   prob sig
## Sign Bias           0.1685 0.8665    
## Negative Sign Bias  0.0495 0.9606    
## Positive Sign Bias  0.5855 0.5594    
## Joint Effect        0.3502 0.9503    
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     33.53     0.020892
## 2    30     52.34     0.004996
## 3    40     54.88     0.047181
## 4    50     74.56     0.010767
## 
## 
## Elapsed time : 1.021805

Q-Q Plot of M1

plot(fit, which=9)

forecast for next 60 months with 99% confidence interval

forecast = ugarchforecast(fit, n.ahead = 60, conf.level = 0.99); forecast
## 
## *------------------------------------*
## *       GARCH Model Forecast         *
## *------------------------------------*
## Model: eGARCH
## Horizon: 60
## Roll Steps: 0
## Out of Sample: 0
## 
## 0-roll forecast [T0=1970-04-29]:
##         Series   Sigma
## T+1  1.034e-02 0.04073
## T+2  4.624e-03 0.03496
## T+3  1.167e-03 0.04102
## T+4  2.944e-04 0.04169
## T+5  7.429e-05 0.04455
## T+6  1.874e-05 0.04619
## T+7  4.730e-06 0.04811
## T+8  1.193e-06 0.04967
## T+9  3.011e-07 0.05116
## T+10 7.599e-08 0.05246
## T+11 1.917e-08 0.05364
## T+12 4.838e-09 0.05469
## T+13 1.221e-09 0.05563
## T+14 3.080e-10 0.05647
## T+15 7.773e-11 0.05721
## T+16 1.961e-11 0.05786
## T+17 4.949e-12 0.05845
## T+18 1.249e-12 0.05896
## T+19 3.151e-13 0.05941
## T+20 7.951e-14 0.05981
## T+21 2.006e-14 0.06016
## T+22 5.062e-15 0.06047
## T+23 1.277e-15 0.06074
## T+24 3.223e-16 0.06098
## T+25 8.133e-17 0.06119
## T+26 2.052e-17 0.06138
## T+27 5.178e-18 0.06154
## T+28 1.307e-18 0.06168
## T+29 3.297e-19 0.06181
## T+30 8.319e-20 0.06192
## T+31 2.099e-20 0.06201
## T+32 5.297e-21 0.06210
## T+33 1.337e-21 0.06217
## T+34 3.372e-22 0.06223
## T+35 8.509e-23 0.06229
## T+36 2.147e-23 0.06234
## T+37 5.418e-24 0.06238
## T+38 1.367e-24 0.06242
## T+39 3.450e-25 0.06246
## T+40 8.704e-26 0.06248
## T+41 2.196e-26 0.06251
## T+42 5.542e-27 0.06253
## T+43 1.398e-27 0.06255
## T+44 3.529e-28 0.06257
## T+45 8.904e-29 0.06258
## T+46 2.247e-29 0.06260
## T+47 5.669e-30 0.06261
## T+48 1.430e-30 0.06262
## T+49 3.609e-31 0.06263
## T+50 9.107e-32 0.06263
## T+51 2.298e-32 0.06264
## T+52 5.799e-33 0.06265
## T+53 1.463e-33 0.06265
## T+54 3.692e-34 0.06266
## T+55 9.316e-35 0.06266
## T+56 2.351e-35 0.06266
## T+57 5.931e-36 0.06267
## T+58 1.497e-36 0.06267
## T+59 3.777e-37 0.06267
## T+60 9.529e-38 0.06267

Plot Forecast

plot(forecast, which = 1, main = "Forecast of Log Returns", xlab = "Time", ylab = "Log Returns")

Write R codes for Kaplan Meier models, results plots usng aml data in the survival package

Load Necessary Libraries

library(survival)
library(ggplot2)
library(survminer)
## Loading required package: ggpubr
## 
## Attaching package: 'survminer'
## The following object is masked from 'package:survival':
## 
##     myeloma
library(ggfortify)
library(survival)

Load the aml dataset

aml_data = aml;aml_data
##    time status             x
## 1     9      1    Maintained
## 2    13      1    Maintained
## 3    13      0    Maintained
## 4    18      1    Maintained
## 5    23      1    Maintained
## 6    28      0    Maintained
## 7    31      1    Maintained
## 8    34      1    Maintained
## 9    45      0    Maintained
## 10   48      1    Maintained
## 11  161      0    Maintained
## 12    5      1 Nonmaintained
## 13    5      1 Nonmaintained
## 14    8      1 Nonmaintained
## 15    8      1 Nonmaintained
## 16   12      1 Nonmaintained
## 17   16      0 Nonmaintained
## 18   23      1 Nonmaintained
## 19   27      1 Nonmaintained
## 20   30      1 Nonmaintained
## 21   33      1 Nonmaintained
## 22   43      1 Nonmaintained
## 23   45      1 Nonmaintained

Fit Kaplan-Meier Survival Model

km_fit = survfit(Surv(time, status) ~ 1, data = aml_data); km_fit
## Call: survfit(formula = Surv(time, status) ~ 1, data = aml_data)
## 
##       n events median 0.95LCL 0.95UCL
## [1,] 23     18     27      18      45

Plot Kaplan-Meier Survival Curve

ggsurvplot(km_fit, data = aml_data, 
           title = "Kaplan-Meier Survival Curve for AML Data",
           xlab = "Time (days)", ylab = "Survival Probability",
           risk.table = TRUE, pval = TRUE)
## Warning in .pvalue(fit, data = data, method = method, pval = pval, pval.coord = pval.coord, : There are no survival curves to be compared. 
##  This is a null model.

Question 4: fit claims data to Weillbull distribution to estimate parameters of distance models between real returns and squared returns for maximum goodness of fit

Load Data & Recquired Libraries

library(fitdistrplus)
claims_data = data$CLAIMS;claims_data
##   [1]  4.59  4.53  4.53  4.63  4.95  5.20  5.61  6.24  6.89  7.40  7.79  7.98
##  [13]  8.15  8.33  8.40  8.36  8.21  8.13  7.84  7.20  6.76  6.43  6.26  6.30
##  [25]  6.43  6.48  6.50  6.47  6.44  6.46  6.59  6.72  6.88  6.87  6.77  6.58
##  [37]  6.42  6.31  6.29  6.34  6.54  6.63  6.65  6.69  6.63  6.63  6.74  6.88
##  [49]  6.97  7.08  7.19  7.33  7.19  7.05  6.85  6.58  6.39  6.21  6.01  5.72
##  [61]  5.39  5.05  4.75  4.50  4.44  4.56  4.96  5.61  6.33  7.24  8.20  9.38
##  [73] 10.67 12.04 13.29 14.33 15.27 15.97 16.40 16.50 16.45 15.93 15.10 14.02
##  [85] 12.82 11.49 10.18  9.00  7.88  6.88  5.96  5.20  4.49  4.05  3.93  3.96
##  [97]  4.02  4.12  4.40  4.69  5.03  5.43  5.85  6.32  7.03  7.88  8.64  9.24
## [109] 10.24 11.42 12.41 13.42 14.35 15.11 15.93 16.72 17.07 16.87 16.56

Fit Weibull Distribution to Claims Data

Kolmogorov-Smirnov (KS) test

fitlnMGEKS <- fitdist(data = claims_data,distr = "weibull", method="mge", gof="KS") 
summary(fitlnMGEKS) 
## Fitting of the distribution ' weibull ' by maximum goodness-of-fit 
## Parameters : 
##       estimate
## shape 2.713481
## scale 8.244151
## Loglikelihood:  -324.4384   AIC:  652.8768   BIC:  658.435

Cramer-von-Mises (CvM)

fitlnMGECvM <- fitdist(data = claims_data,distr = "weibull", method="mge", gof="CvM") 
summary(fitlnMGECvM)
## Fitting of the distribution ' weibull ' by maximum goodness-of-fit 
## Parameters : 
##       estimate
## shape  3.99628
## scale  7.71185
## Loglikelihood:  -488.2004   AIC:  980.4008   BIC:  985.9591

Anderson-Darling (AD)

fitlnMGEAD <- fitdist(data = claims_data,distr = "weibull", method="mge", gof="AD") 
summary(fitlnMGEAD)
## Fitting of the distribution ' weibull ' by maximum goodness-of-fit 
## Parameters : 
##       estimate
## shape 2.260180
## scale 8.915078
## Loglikelihood:  -313.8128   AIC:  631.6256   BIC:  637.1839