Search This Blog

Friday, March 16, 2012

R PROGRAMMING: part#2

Consider 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) In previous assignment you have observed the superimposition of the normal density
separately for SP-rtn and C-rtn data. Now use your intuition or play trial and error method
1
to fit the same data sets with some t-distribution, double-exponential, cauchy as well as a
mixture of two normal distributions on the same graph. Which one is coming best ? What
is you intuition while making the trial and error for mixture distribution ? What is your
comment on the tail of cauchy distribution (something odd !!!! pay more attention !!!) (You
may illustrate the answer after making part (b) and part (c) too).
(b) Make different qqplots taking quantiles for each density you have chosen against the
empirical quantiles of the data. Give proper interpretation.
(c) Plot the empirical survival curve and superimpose survival curves for all the densities
you have chosen to fit the data. Comment on the picture.
Extra : If you find the task really interesting; Play a little more with mixture of two
cauchy distributions too.


Conclusion

Which one is coming best ?
ANS:- As I see Normal distribution has havier tail than all others and Cauchy Distribution has lighter tail than all others. Similar for their mixture distribution but tail of mixture distribution is slightly lighter than actual normal distribution and similarly tail of mixture of cauchy distribution is slightly lighter than actual cauchy distribution.
According to my observation Mixture of two normal distribution is best fitted distribution with our C_rtn and SP_rtn data, hence it is best.

What is you intuition while making the trial and error for mixture distribution ?
ANS:- I have taken 90% of N(0,1) and 10% of N(0,25) for mixture of two normal distribution and 90% of C(0,1) and 10% of C(0,5) for mixture of two Cauchy distribution. Variance of mixture normal dis is 3.4, we get different distribution of mixture of these two and N(0,3.4) although both distribution have same mean and variance. Same for Cauchy Distribution.

What is your comment on the tail of cauchy distribution (something odd !!!! pay more attention !!!)
ANS:- Tail of Cauchy distribution is lighter than all other distribution drown and Tail of mixture of two cauchy distribution is slightly lighter than actual cauchy distribution.

Q a >>>>> Code for C_rtn

crtn=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];
mean_C=mean(C);
mean_SP=mean(SP);
var_C=var(C);
var_SP=var(SP);
c_rtn=t(C);
sp_rtn=t(SP);
# 4 figures arranged in 2 rows and 2 columns
#attach(mtcars)
par(mfrow=c(2,3))
###---------- Normal Plot- ------------####
h1=hist(c_rtn, breaks=400, col="red", xlab="c_rtn", freq=FALSE, xlim=c(-0.2,0.2), ylim=c(0,40))
x=seq(min(c_rtn),max(c_rtn),length=200)
ynorm=dnorm(x, mean=0, sd=0.012);
#ynorm=ynorm*diff(h1$mids[1:2])*length(c_rtn);
par(new=TRUE)
plot(x, ynorm, type="l", col="blue", lwd=2 , xlim=c(-0.2,0.2), ylim=c(0,40));
labels=c("normal Distribution ")
legend("topleft", border="white", inset=0, title="Distributions", labels, lwd=2, lty=c(1,1,1,1,2))
####--------student T plot ------------#######
h1=hist(c_rtn, breaks=400, col="red", xlab="c_rtn", freq=FALSE, xlim=c(-0.2,0.2), ylim=c(0,40))
x=seq(-15,15,length=200);
x1=seq(-0.2,0.2,length=200);
yt=dt(x, 200000);
yt=yt*diff(h1$mids[1:2])*length(c_rtn)*20;
par(new=TRUE)
plot(x1,yt, type="l", col="blue", lwd=2, axes=FALSE, xlim=c(-0.2,0.2), ylim=c(0,40))
labels=c("Student t Distribution ")
legend("topleft", border="white", inset=0, title="Distributions", labels, lwd=2, lty=c(1,1,1,1,2))
#####---------------- Double exponential or Laplace plot ----------##########
h1=hist(c_rtn, breaks=400, col="red", xlab="c_rtn", freq=FALSE, xlim=c(-0.2,0.2), ylim=c(0,40))
b=1;
mu=0;
ylap=(1/(2*b))*(exp(-abs(x-mu)/b))
ylap=ylap*diff(h1$mids[1:2])*length(c_rtn)*18;
par(new=TRUE)
plot(x1,ylap, type="l", col="blue", lwd=2, axes=FALSE, xlim=c(-0.2,0.2), ylim=c(0,40))
labels=c("laplace Distribution ")
legend("topleft", border="white", inset=0, title="Distributions", labels, lwd=2, lty=c(1,1,1,1,2))
######--------------- Cauchy Distribution ---------- #########
h1=hist(c_rtn, breaks=400, col="red", xlab="c_rtn", freq=FALSE, xlim=c(-0.2,0.2), ylim=c(0,40))
x=seq(-15,15,length=200);
x1=seq(-0.2,0.2,length=200);
gama=0.75;
x0=0;
yc=dcauchy(x,location=x0, scale=gama, log=FALSE)
yc=yc*diff(h1$mids[1:2])*length(c_rtn)*20;
par(new=TRUE)
plot(x1,yc, type="l", col="blue", lwd=2, axes=FALSE, xlim=c(-0.2,0.2), ylim=c(0,40))
labels=c("Cauchy Distribution ")
legend("topleft", border="white", inset=0, title="Distributions", labels, lwd=2, lty=c(1,1,1,1,2))
########---------- Mixture of two normal distribution -------------########
h1=hist(c_rtn, breaks=400, col="red", xlab="c_rtn", freq=FALSE, xlim=c(-0.2,0.2), ylim=c(0,40))
x=seq(min(c_rtn),max(c_rtn),length=200)
ymnorm=0.9*dnorm(x,mean=0, sd=0.012) + 0.1*dnorm(x,mean=0, sd=0.009)
#ymnorm=ymnorm*diff(h1$mids[1:2])*length(c_rtn)*20;
par(new=TRUE)
plot(x,ymnorm, type="l", col="blue", lwd=2, axes=FALSE, xlim=c(-0.2,0.2), ylim=c(0,40))
labels=c("mixture of normal Distribution ")
legend("topleft", border="white", inset=0, title="Distributions", labels, lwd=2, lty=c(1,1,1,1,2))
######------------- Mixture of two Cauchy Distribution ------------- ########
h1=hist(c_rtn, breaks=400, col="red", xlab="c_rtn", freq=FALSE, xlim=c(-0.2,0.2), ylim=c(0,40))
x=seq(-15,15,length=200);
x1=seq(-0.2,0.2,length=200);
ymc=0.9*dcauchy(x,location=0, scale=0.75, log=FALSE) + 0.1*dcauchy(x,location=0, scale=0.9, log=FALSE)
ymc=ymc*diff(h1$mids[1:2])*length(c_rtn)*20;
par(new=TRUE)
plot(x1,ymc, type="l", col="blue", lwd=2, axes=FALSE, xlim=c(-0.2,0.2), ylim=c(0,40))
labels=c("Mixture of cauchy Distribution ")
legend("topleft", border="white", inset=0, title="Distributions", labels, lwd=2, lty=c(1,1,1,1,2))
########## Labels and legends ###################
#labels=c("normal = points", "student t = lines", "laplace = both", "cauchy = both overplotted", "Mixture of normal = steps", "Mixture of cauchy = other steps")
#legend("topleft", border="white", inset=0, title="Distributions", labels, lwd=2, lty=c(1,1,1,1,2))
#title(" Distribution Plots of SP_rtn with different density functions")
}


Q a >>>>> Code for SP_rtn

sprtn=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];
mean_C=mean(C);
mean_SP=mean(SP);
var_C=var(C);
var_SP=var(SP);
c_rtn=t(C);
sp_rtn=t(SP);
# 4 figures arranged in 2 rows and 2 columns
#attach(mtcars)
par(mfrow=c(2,3))
###---------- Normal Plot- ------------####
h1=hist(sp_rtn, breaks=400, col="red", xlab="sp_rtn", freq=FALSE, xlim=c(-0.1,0.1), ylim=c(0,80))
x=seq(min(sp_rtn),max(sp_rtn),length=200)
ynorm=dnorm(x, mean=0, sd=0.0065);
#ynorm=ynorm*diff(h1$mids[1:2])*length(c_rtn);
par(new=TRUE)
plot(x, ynorm, type="l", col="blue", lwd=2 , xlim=c(-0.1,0.1), ylim=c(0,80));
labels=c("normal Distribution ")
legend("topleft", border="white", inset=0, title="Distributions", labels, lwd=2, lty=c(1,1,1,1,2))
####--------student T plot ------------#######
h1=hist(sp_rtn, breaks=400, col="red", xlab="sp_rtn", freq=FALSE, xlim=c(-0.1,0.1), ylim=c(0,80))
x=seq(-15,15,length=200);
x1=seq(-0.1,0.1,length=200);
yt=dt(x, 200000);
yt=yt*diff(h1$mids[1:2])*length(c_rtn)*160;
par(new=TRUE)
plot(x1,yt, type="l", col="blue", lwd=2, axes=FALSE, xlim=c(-0.1,0.1), ylim=c(0,80))
labels=c("Student t Distribution ")
legend("topleft", border="white", inset=0, title="Distributions", labels, lwd=2, lty=c(1,1,1,1,2))
#####---------------- Double exponential or Laplace plot ----------##########
h1=hist(sp_rtn, breaks=400, col="red", xlab="sp_rtn", freq=FALSE, xlim=c(-0.1,0.1), ylim=c(0,80))
b=1;
mu=0;
ylap=(1/(2*b))*(exp(-abs(x-mu)/b))
ylap=ylap*diff(h1$mids[1:2])*length(c_rtn)*135;
par(new=TRUE)
plot(x1,ylap, type="l", col="blue", lwd=2, axes=FALSE, xlim=c(-0.1,0.1), ylim=c(0,80))
labels=c("laplace Distribution ")
legend("topleft", border="white", inset=0, title="Distributions", labels, lwd=2, lty=c(1,1,1,1,2))
######--------------- Cauchy Distribution ---------- #########
h1=hist(sp_rtn, breaks=400, col="red", xlab="sp_rtn", freq=FALSE, xlim=c(-0.1,0.1), ylim=c(0,80))
x=seq(-15,15,length=200);
x1=seq(-0.1,0.1,length=200);
gama=0.75;
x0=0;
yc=dcauchy(x,location=x0, scale=gama, log=FALSE)
yc=yc*diff(h1$mids[1:2])*length(c_rtn)*150;
par(new=TRUE)
plot(x1,yc, type="l", col="blue", lwd=2, axes=FALSE, xlim=c(-0.1,0.1), ylim=c(0,80))
labels=c("Cauchy Distribution ")
legend("topleft", border="white", inset=0, title="Distributions", labels, lwd=2, lty=c(1,1,1,1,2))
########---------- Mixture of two normal distribution -------------########
h1=hist(sp_rtn, breaks=400, col="red", xlab="sp_rtn", freq=FALSE, xlim=c(-0.1,0.1), ylim=c(0,80))
x=seq(min(sp_rtn),max(sp_rtn),length=200)
ymnorm=0.9*dnorm(x,mean=0, sd=0.0061) + 0.1*dnorm(x,mean=0, sd=0.0099)
#ymnorm=ymnorm*diff(h1$mids[1:2])*length(c_rtn)*20;
par(new=TRUE)
plot(x,ymnorm, type="l", col="blue", lwd=2, axes=FALSE, xlim=c(-0.1,0.1), ylim=c(0,80))
labels=c("mixture of normal Distribution ")
legend("topleft", border="white", inset=0, title="Distributions", labels, lwd=2, lty=c(1,1,1,1,2))
######------------- Mixture of two Cauchy Distribution ------------- ########
h1=hist(sp_rtn, breaks=400, col="red", xlab="sp_rtn", freq=FALSE, xlim=c(-0.1,0.1), ylim=c(0,80))
x=seq(-15,15,length=200);
x1=seq(-0.1,0.1,length=200);
ymc=0.9*dcauchy(x,location=0, scale=0.75, log=FALSE) + 0.1*dcauchy(x,location=0, scale=0.9, log=FALSE)
ymc=ymc*diff(h1$mids[1:2])*length(c_rtn)*150;
par(new=TRUE)
plot(x1,ymc, type="l", col="blue", lwd=2, axes=FALSE, xlim=c(-0.1,0.1), ylim=c(0,80))
labels=c("Mixture of cauchy Distribution ")
legend("topleft", border="white", inset=0, title="Distributions", labels, lwd=2, lty=c(1,1,1,1,2))
########## Labels and legends ###################
#labels=c("normal = points", "student t = lines", "laplace = both", "cauchy = both overplotted", "Mixture of normal = steps", "Mixture of cauchy = other steps")
#legend("topleft", border="white", inset=0, title="Distributions", labels, lwd=2, lty=c(1,1,1,1,2))
#title(" Distribution Plots of SP_rtn with different density functions")
}


Q b >>>>>>>>>>>> Code for C_rtn

qbc=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];
mean_C=mean(C);
mean_SP=mean(SP);
var_C=var(C);
var_SP=var(SP);
c_rtn=t(C);
sp_rtn=t(SP);

# 4 figures arranged in 2 rows and 2 columns
#attach(mtcars)
par(mfrow=c(2,3))
x=seq(-10,10,length=10000);
qc_rtn=quantile(c_rtn,prob=seq(0,1,.0001))
###---------- Normal Plot- ------------####
yqnorm=qnorm(x, mean=0, sd=1, lower.tail = TRUE, log.p = FALSE);
qqplot(yqnorm,qc_rtn)
qqline(qc_rtn)
labels=c("Normal Distribution")
legend("topright", inset=.05, title="Q-Q Plots c_rtn", labels, lwd=2)
####--------student T plot ------------#######
yqt=qt(x, 3, lower.tail = TRUE, log.p = FALSE);
qqplot(yqt,qc_rtn)
qqline(qc_rtn)
labels=c("Student t Distribution")
legend("topright", inset=.05, title="Q-Q Plots c_rtn ", labels, lwd=2)
#####---------------- Double exponential or Laplace plot ----------##########
#ylap=dlaplace(x,location=0,scale=1)
b=1;
mu=0;
ylap=(1/(2*b))*(exp(-abs(x-mu)/b))
yqlap=quantile(ylap)
qqplot(yqlap,qc_rtn)
qqline(qc_rtn)
labels=c("Laplace Distribution")
legend("topright", inset=.05, title="Q-Q Plots c_rtn ", labels, lwd=2)
######--------------- Cauchy Distribution ---------- #########
yqcauchy=qcauchy(x, location=0, scale=1, lower.tail = TRUE, log.p = FALSE);
qqplot(yqcauchy,qc_rtn)
qqline(qc_rtn)
labels=c("Cauchy Distribution")
legend("topright", inset=.05, title="Q-Q Plots c_rtn ", labels, lwd=2)
########---------- Mixture of two normal distribution -------------########
yqmnorm=0.9*qnorm(x,mean=0, sd=1) + 0.1*qnorm(x,mean=0, sd=5)
qqplot(yqmnorm,qc_rtn)
qqline(qc_rtn)
labels=c("Mixture of Two Normal")
legend("topright", inset=.05, title="Q-Q Plots c_rtn ", labels, lwd=2)
######------------- Mixture of two Cauchy Distribution ------------- ########
yqmc=0.9*qcauchy(x,location=0, scale=1, log=FALSE) + 0.1*qcauchy(x,location=0, scale=5, log=FALSE)
qqplot(yqmc,qc_rtn)
qqline(qc_rtn)
labels=c("Mixture of Two cauchy")
legend("topright", inset=.05, title="Q-Q Plots c_rtn ", labels, lwd=2)
########## Labels and legends ###################
#labels=c("normal = p", "student t = l", "laplace = b", "cauchy = o", "Mixture of normal = s", "Mixture of cauchy = S")
#legend("topright", inset=.05, title="Distributions", labels, lwd=2, lty=c(1,1,1,1,2))
}



Q b >>>>> code for SP_rtn

qbsp=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];
mean_C=mean(C);
mean_SP=mean(SP);
var_C=var(C);
var_SP=var(SP);
c_rtn=t(C);
sp_rtn=t(SP);

# 4 figures arranged in 2 rows and 2 columns
#attach(mtcars)
par(mfrow=c(2,3))
x=seq(-10,10,length=10000);
qsp_rtn=quantile(sp_rtn,prob=seq(0,1,.0001))
###---------- Normal Plot- ------------####
yqnorm=qnorm(x, mean=0, sd=1, lower.tail = TRUE, log.p = FALSE);
qqplot(yqnorm,qsp_rtn)
qqline(qsp_rtn)
labels=c("Normal Distribution")
legend("topright", inset=.05, title="Q-Q Plots sp_rtn", labels, lwd=2)
####--------student T plot ------------#######
yqt=qt(x, 3, lower.tail = TRUE, log.p = FALSE);
qqplot(yqt,qsp_rtn)
qqline(qsp_rtn)
labels=c("Student t Distribution")
legend("topright", inset=.05, title="Q-Q Plots sp_rtn", labels, lwd=2)
#####---------------- Double exponential or Laplace plot ----------##########
#ylap=dlaplace(x,location=0,scale=1)
b=1;
mu=0;
ylap=(1/(2*b))*(exp(-abs(x-mu)/b))
yqlap=quantile(ylap)
qqplot(yqlap,qsp_rtn)
qqline(qsp_rtn)
labels=c("Laplace Distribution")
legend("topright", inset=.05, title="Q-Q Plots sp_rtn", labels, lwd=2)
######--------------- Cauchy Distribution ---------- #########
yqcauchy=qcauchy(x, location=0, scale=1, lower.tail = TRUE, log.p = FALSE);
qqplot(yqcauchy,qsp_rtn)
qqline(qsp_rtn)
labels=c("Cauchy Distribution")
legend("topright", inset=.05, title="Q-Q Plots sp_rtn", labels, lwd=2)
########---------- Mixture of two normal distribution -------------########
yqmnorm=0.9*qnorm(x,mean=0, sd=1) + 0.1*qnorm(x,mean=0, sd=5)
qqplot(yqmnorm,qsp_rtn)
qqline(qsp_rtn)
labels=c("Mixture of Two Normal")
legend("topright", inset=.05, title="Q-Q Plots sp_rtn", labels, lwd=2)
######------------- Mixture of two Cauchy Distribution ------------- ########
yqmc=0.9*qcauchy(x,location=0, scale=1, log=FALSE) + 0.1*qcauchy(x,location=0, scale=5, log=FALSE)
qqplot(yqmc,qsp_rtn)
qqline(qsp_rtn)
labels=c("Mixture of Two cauchy")
legend("topright", inset=.05, title="Q-Q Plots sp_rtn", labels, lwd=2)
########## Labels and legends ###################
#labels=c("normal = p", "student t = l", "laplace = b", "cauchy = o", "Mixture of normal = s", "Mixture of cauchy = S")
#legend("topright", inset=.05, title="Distributions", labels, lwd=2, lty=c(1,1,1,1,2))
}



Q c >>>>>>>> Code for C_rtn

qcc=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];
mean_C=mean(C);
mean_SP=mean(SP);
var_C=var(C);
var_SP=var(SP);
c_rtn=t(C);
sp_rtn=t(SP);

# 4 figures arranged in 2 rows and 2 columns
#attach(mtcars)
par(mfrow=c(3,3))
#x=seq(-10,10,length=1000);
#x=seq(min(c_rtn), max(c_rtn), length=length(c_rtn))
h1=hist(c_rtn, breaks=400, col="red", xlab="c_rtn", freq=FALSE, xlim=c(-0.2,0.2), ylim=c(0,40))
###---------- Normal Plot- ------------####
x=seq(min(c_rtn),max(c_rtn),length=200)
ydnorm=dnorm(x, mean=0, sd=0.012);
plot.stepfun(ydnorm, col.vert = "gray20", main = "empirical survival curve c_rtn of Normal distribution ", axes=FALSE, xlim=c(-.2,.6) , col="red")
par(new=TRUE)
plot.stepfun(c_rtn,col.vert="gray20" , main = "empirical survival curve c_rtn of Normal distribution " , xlim=c(-.2,.6), col="blue");

###--------student T plot ------------#######
x=seq(-15,15,length=200);
x1=seq(-0.2,0.2,length=200);
yt=dt(x, 200000);
ydt=yt*diff(h1$mids[1:2])*length(c_rtn)*20;
plot.stepfun(ydt, col.vert = "gray20", main = "empirical survival curve c_rtn of Student t distribution " , axes=FALSE, xlim=c(-.2,.6), col="red")
par(new=TRUE)
plot.stepfun(c_rtn,col.vert="gray20", main = "empirical survival curve c_rtn of Student i distribution " , xlim=c(-.2,.6) , col="blue");
#####---------------- Double exponential or Laplace plot ----------##########
#ylap=dlaplace(x,location=0,scale=1)
b=1;
mu=0;
ylap=(1/(2*b))*(exp(-abs(x-mu)/b))
ylap=ylap*diff(h1$mids[1:2])*length(c_rtn)*18;
plot.stepfun(ylap, col.vert = "gray20", main = "empirical survival curve c_rtn of Laplace distribution " , axes=FALSE, xlim=c(-.2,.6) , col="red")
par(new=TRUE)
plot.stepfun(c_rtn,col.vert="gray20", main = "empirical survival curve c_rtn of Laplace distribution " , xlim=c(-.2,.6) , col="blue" );

######--------------- Cauchy Distribution ---------- #########
x=seq(-15,15,length=200);
x1=seq(-0.2,0.2,length=200);
gama=0.75;
x0=0;
yc=dcauchy(x,location=x0, scale=gama, log=FALSE)
ydc=yc*diff(h1$mids[1:2])*length(c_rtn)*20;
plot.stepfun(ydc, col.vert = "gray20", main = "empirical survival curve c_rtn of Cauchy distribution " , axes=FALSE, xlim=c(-.2,.6) , col="red")
par(new=TRUE)
plot.stepfun(c_rtn,col.vert="gray20", main = "empirical survival curve c_rtn of Cauchy distribution " , xlim=c(-.2,.6) , col="blue");
########---------- Mixture of two normal distribution -------------########
x=seq(min(c_rtn),max(c_rtn),length=200)
ydmnorm=0.9*dnorm(x,mean=0, sd=0.012) + 0.1*dnorm(x,mean=0, sd=0.009)
plot.stepfun(ydmnorm, col.vert = "gray20", main = "empirical survival curve c_rtn of Mixture of normal distribution " , axes=FALSE, xlim=c(-.2,.6) , col="red")
par(new=TRUE)
plot.stepfun(c_rtn,col.vert="gray20", main = "empirical survival curve c_rtn of Mixture of Normal distribution " , xlim=c(-.2,.6) , col="blue");
######------------- Mixture of two Cauchy Distribution ------------- ########
x=seq(-15,15,length=200);
x1=seq(-0.2,0.2,length=200);
ymc=0.9*dcauchy(x,location=0, scale=0.75, log=FALSE) + 0.1*dcauchy(x,location=0, scale=0.9, log=FALSE)
ydmc=ymc*diff(h1$mids[1:2])*length(c_rtn)*20;
plot.stepfun(ydmc, col.vert = "gray20", main = "empirical survival curve c_rtn of Mixture of cauchy distribution " , axes=FALSE, xlim=c(-.2,.6) , col="red")
par(new=TRUE)
plot.stepfun(c_rtn,col.vert="gray20", main = "empirical survival curve c_rtn of Mixture of cauchy distribution " , xlim=c(-.2,.6) , col="blue");
########## Labels and legends ###################
#labels=c("normal = p", "student t = l", "laplace = b", "cauchy = o", "Mixture of normal = s", "Mixture of cauchy = S")
#legend("topright", inset=.05, title="Distributions", labels, lwd=2, lty=c(1,1,1,1,2))
}




Q c >>>>>>>> Code for SP_rtn

qcsp=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];
mean_C=mean(C);
mean_SP=mean(SP);
var_C=var(C);
var_SP=var(SP);
c_rtn=t(C);
sp_rtn=t(SP);

# 4 figures arranged in 2 rows and 2 columns
#attach(mtcars)
par(mfrow=c(3,3))
#x=seq(-10,10,length=1000);
#x=seq(min(sp_rtn), max(sp_rtn), length=length(sp_rtn))
h1=hist(sp_rtn, breaks=400, col="red", xlab="sp_rtn", freq=FALSE, xlim=c(-0.1,0.1), ylim=c(0,80))
###---------- Normal Plot- ------------####
x=seq(min(sp_rtn),max(sp_rtn),length=200)
ydnorm=dnorm(x, mean=0, sd=0.0065);
plot.stepfun(ydnorm, col.vert = "gray20", main = "empirical survival curve SP_rtn of Normal distribution ", axes=FALSE, xlim=c(-.2,.6) , col="red")
par(new=TRUE)
plot.stepfun(sp_rtn,col.vert="gray20" , main = "empirical survival curve SP_rtn of Normal distribution " , xlim=c(-.2,.6), col="blue");

###--------student T plot ------------#######
x=seq(-15,15,length=200);
x1=seq(-0.1,0.1,length=200);
yt=dt(x, 200000);
ydt=yt*diff(h1$mids[1:2])*length(c_rtn)*160;
plot.stepfun(ydt, col.vert = "gray20", main = "empirical survival curve sp_rtn of Student t distribution " , axes=FALSE, xlim=c(-.2,.6), col="red")
par(new=TRUE)
plot.stepfun(sp_rtn,col.vert="gray20", main = "empirical survival curve sp_rtn of Student i distribution " , xlim=c(-.2,.6) , col="blue");
#####---------------- Double exponential or Laplace plot ----------##########
#ylap=dlaplace(x,location=0,scale=1)
b=1;
mu=0;
ylap=(1/(2*b))*(exp(-abs(x-mu)/b))
ylap=ylap*diff(h1$mids[1:2])*length(c_rtn)*135;
plot.stepfun(ylap, col.vert = "gray20", main = "empirical survival curve sp_rtn of Laplace distribution " , axes=FALSE, xlim=c(-.2,.6) , col="red")
par(new=TRUE)
plot.stepfun(sp_rtn,col.vert="gray20", main = "empirical survival curve sp_rtn of Laplace distribution " , xlim=c(-.2,.6) , col="blue" );

######--------------- Cauchy Distribution ---------- #########
x=seq(-15,15,length=200);
x1=seq(-0.1,0.1,length=200);
gama=0.75;
x0=0;
yc=dcauchy(x,location=x0, scale=gama, log=FALSE)
ydc=yc*diff(h1$mids[1:2])*length(c_rtn)*150;
plot.stepfun(ydc, col.vert = "gray20", main = "empirical survival curve sp_rtn of Cauchy distribution " , axes=FALSE, xlim=c(-.2,.6) , col="red")
par(new=TRUE)
plot.stepfun(sp_rtn,col.vert="gray20", main = "empirical survival curve sp_rtn of Cauchy distribution " , xlim=c(-.2,.6) , col="blue");
########---------- Mixture of two normal distribution -------------########
x=seq(min(sp_rtn),max(sp_rtn),length=200)
ydmnorm=0.9*dnorm(x,mean=0, sd=0.0061) + 0.1*dnorm(x,mean=0, sd=0.0099)
plot.stepfun(ydmnorm, col.vert = "gray20", main = "empirical survival curve sp_rtn of Mixture of normal distribution " , axes=FALSE, xlim=c(-.2,.6) , col="red")
par(new=TRUE)
plot.stepfun(sp_rtn,col.vert="gray20", main = "empirical survival curve sp_rtn of Mixture of Normal distribution " , xlim=c(-.2,.6) , col="blue");
######------------- Mixture of two Cauchy Distribution ------------- ########
x=seq(-15,15,length=200);
x1=seq(-0.1,0.1,length=200);
ymc=0.9*dcauchy(x,location=0, scale=0.75, log=FALSE) + 0.1*dcauchy(x,location=0, scale=0.9, log=FALSE)
ydmc=ymc*diff(h1$mids[1:2])*length(c_rtn)*150;
plot.stepfun(ydmc, col.vert = "gray20", main = "empirical survival curve sp_rtn of Mixture of cauchy distribution " , axes=FALSE, xlim=c(-.2,.6) , col="red")
par(new=TRUE)
plot.stepfun(sp_rtn,col.vert="gray20", main = "empirical survival curve sp_rtn of Mixture of cauchy distribution " , xlim=c(-.2,.6) , col="blue");
########## Labels and legends ###################
#labels=c("normal = p", "student t = l", "laplace = b", "cauchy = o", "Mixture of normal = s", "Mixture of cauchy = S")
#legend("topright", inset=.05, title="Distributions", labels, lwd=2, lty=c(1,1,1,1,2))
}



********* End ********

No comments:

Post a Comment

Thank you