R Programming: part#3

onsider the daily stock return of the Citigroup (tick symbol C) and the Standard and Poor’s
500 Composite index from January 2001 to December 2008. The data are simple returns
and in the file d-csp0108.txt (three columns with date, C-rtn, SP-rtn).
(a) Make kernel density plot for two data sets based on second order epanechnikov and
gaussian kernel. Use Silverman’s thumb rule for the choice of bandwidth. Comment on your
observation.
1
(b) Transform simple returns into log-returns. Make the kernel density plot based on
second order epanechnikov and gaussian kernel. Also fit log-normal distribution taking
parameters as maximum likelihood estimates.
(c) Draw a random sample of size 20, 40, 100, 200 from Weibull(shape para=2,scale
para=1). Calculate maximum likelihood estimates based on all samples. What is your
observation ?

R/S-Plus CODE

Q No 3a>>>>>>>>>>>>>>>>

kernel=function()

{

####### --------- Lybraries ------------#####

library("MASS");

require(graphics);

data_set=read.table("d-csp0108.txt",header=TRUE);

date=data_set[1];

C=data_set[2];

SP=data_set[3];

c_rtn=t(C);

sp_rtn=t(SP);

# 4 figures arranged in 2 rows and 2 columns

#attach(mtcars)

par(mfrow=c(2,2))

######### KDE plot for c_rtn For #############

cgd=density(c_rtn, bw.nrd0(c(c_rtn)), adjust=1, kernel=c("gaussian") );

ced=density(c_rtn, bw.nrd0(c(c_rtn)), adjust=1, kernel=c("epanechnikov") );

plot(cgd, type="l", col="red", lwd = 2)

par(new=TRUE)

plot(ced, type="p", col="blue", lwd=2)

labels=c("KDE Plot of Gaussian kernel", "KDE of epanechnikov kernel ")

legend("topright", border="white", inset=0, title="KDE of C_rtn", labels, lwd=2, lty=c(1,1,1,1,2))

plot(cgd, type="s", col="red", lwd = 2, xlim=c(0.01, 0.01001))

par(new=TRUE)

plot(ced, type="S", col="blue", lwd=2, xlim=c(0.01, 0.01001))

labels=c("KDE Zoomed plots of Gaussian Kernel", "KDE Zoomed plots of epanechnikov Kernel")

legend("topright", border="white", inset=0, title="KDE of C_rtn", labels, lwd=2, lty=c(1,1,1,1,2))

######### KDE plot for sp_rtn For #############

spgd=density(sp_rtn, bw.nrd0(c(sp_rtn)), adjust=1, kernel=c("gaussian") );

sped=density(sp_rtn, bw.nrd0(c(sp_rtn)), adjust=1, kernel=c("epanechnikov") );

plot(spgd, type="l", col="red", lwd = 2)

par(new=TRUE)

plot(sped, type="p", col="blue", lwd=2)

labels=c("KDE Plot of Gaussian kernel", "KDE of epanechnikov kernel ")

legend("topright", border="white", inset=0, title="KDE of SP_rtn", labels, lwd=2, lty=c(1,1,1,1,2))

plot(spgd, type="s", col="red", lwd = 2, xlim=c(0.01, 0.01001))

par(new=TRUE)

plot(sped, type="S", col="blue", lwd=2, xlim=c(0.01, 0.01001))

labels=c("KDE Zoomed plots of Gaussian Kernel", "KDE Zoomed plots of epanechnikov Kernel")

legend("topright", border="white", inset=0, title="SP_rtn", labels, lwd=2, lty=c(1,1,1,1,2))

amiseg=(5/4) * ( ( 1.06 * 1 * ( ( 1 / ( 2 * ( sqrt(3.14) ) ) )^4 ) )^(1/5) ) * (2011^(-4/5))

print(amiseg)

amisee=(5/4) * ( ( 2.34 * ( (1/5)^2) * ( ( 3 / 5 )^4 ) )^(1/5) ) * (2011^(-4/5))

print(amisee)

a= ( ( 1 / ( 2 * ( sqrt(3.14) ) ) )^4 ) * 1.06

b=( ( ( 3 / 5 )^4 ) *( (1/5)^2) ) * 2.34

print(a)

print(b)

}

Q No 3b>>>>>>>>>>>>>>>>

logrtn=function()

{

####### --------- Lybraries ------------#####

library("MASS");

require(graphics);

data_set=read.table("d-csp0108.txt",header=TRUE);

date=data_set[1];

C=data_set[2];

SP=data_set[3];

C=log(1+C);

SP=log(1+SP);

c_rtn=t(C);

sp_rtn=t(SP);

# 4 figures arranged in 2 rows and 2 columns

#attach(mtcars)

par(mfrow=c(2,2))

######### KDE plot for c_rtn For #############

cgd=density(c_rtn, bw.nrd0(c(c_rtn)), adjust=1, kernel=c("gaussian") );

ced=density(c_rtn, bw.nrd0(c(c_rtn)), adjust=1, kernel=c("epanechnikov") );

x=seq(min(c_rtn),max(c_rtn),length=200)

ynorm=dnorm(x, mean=0, sd=0.042);

#d=dnorm(x, mean=0, sd=1)

l=exp(ynorm)

fitdistr(l, "lognormal");

x1=seq(-.31455, .46426, length=200)

plot(cgd, type="l", col="red", lwd = 2)

par(new=TRUE)

plot(ced, type="p", col="blue", lwd=2)

par(new=TRUE)

plot(x1, l, type="S", col="green", lwd=2 , xlim=c(-.31455, .46426) );

labels=c("KDE Plot of Gaussian kernel", "KDE of epanechnikov kernel ", "fit lnorm by MLE")

legend("topright", border="white", inset=0, title="log rtn KDE of C_rtn", labels, lwd=2, lty=c(1,1,1,1,2))

plot(cgd, type="s", col="red", lwd = 2, xlim=c(0.01, 0.01001))

par(new=TRUE)

plot(ced, type="S", col="blue", lwd=2, xlim=c(0.01, 0.01001))

labels=c("KDE Zoomed plots of Gaussian Kernel", "KDE Zoomed plots of epanechnikov Kernel")

legend("topright", border="white", inset=0, title="log rtn KDE of C_rtn", labels, lwd=2, lty=c(1,1,1,1,2))

######### KDE plot for sp_rtn For #############

spgd=density(sp_rtn, bw.nrd0(c(sp_rtn)), adjust=1, kernel=c("gaussian") );

sped=density(sp_rtn, bw.nrd0(c(sp_rtn)), adjust=1, kernel=c("epanechnikov") );

xs=seq(min(sp_rtn),max(sp_rtn),length=200)

ynorms=dnorm(xs, mean=0, sd=0.03);

#d=dnorm(x, mean=0, sd=1)

ls=exp(ynorms)

fitdistr(ls, "lognormal");

x1s=seq(-.099964, .11484, length=200)

plot(spgd, type="l", col="red", lwd = 2)

par(new=TRUE)

plot(sped, type="p", col="blue", lwd=2)

par(new=TRUE)

plot(x1s, ls, type="S", col="green", lwd=2 , xlim=c(-.099964, .11484) );

labels=c("KDE Plot of Gaussian kernel", "KDE of epanechnikov kernel ", "fit lnorm by MLE")

legend("topright", border="white", inset=0, title="log rtn KDE of SP_rtn", labels, lwd=2, lty=c(1,1,1,1,2))

plot(spgd, type="s", col="red", lwd = 2, xlim=c(0.01, 0.01001))

par(new=TRUE)

plot(sped, type="S", col="blue", lwd=2, xlim=c(0.01, 0.01001))

labels=c("KDE Zoomed plots of Gaussian Kernel", "KDE Zoomed plots of epanechnikov Kernel")

legend("topright", border="white", inset=0, title="log rtn KDE of SP_rtn", labels, lwd=2, lty=c(1,1,1,1,2))

}

Q No 3c>>>>>>>>>>>>>>>>>

weibull=function(n)

{

####### --------- Lybraries ------------#####

library("MASS");

require(graphics);

x=0;

y=0;

shape=0;

scale=0

a=0;

b=0;

for(i in 1:n)

{

x=rweibull(i+1, shape = 2, scale = 1)

y=fitdistr(x, "weibull")

shape[i]=y$estimate["shape"];

scale[i]=y$estimate["scale"];

}

for(i in 1:n)

{

a[i]=(shape[i]-2)^2;

b[i]=(scale[i]-1)^2;

}

mse_shape=mean(a);

mse_scale=mean(b);

print(mse_shape);

print(mse_scale);

}

***************** END ***************