2014년 4월 12일 토요일

마르코브

library(stringr)

setwd("C:/rdata")
source("DTMC.R")

###here kospi data
dat <- read.csv('port3.csv', header=TRUE, na="")

odat<-dat[c(1:nrow(dat)),]

attach(odat)

x<-odat[order(일자),]


price <- as.numeric(str_replace_all(x$리홈쿠첸, ",", ""))


d <- diff(price)

up2up <-0
up2down <-0
down2up <-0
down2down <-0

for (i in 2:length(d)){
    if (d[i-1]> 0 && d[i]>0){
    up2up <- up2up+1
    }
    if (d[i-1]> 0 && d[i]<0){
    up2down <- up2down +1
    }
    if (d[i-1]< 0 && d[i]>0){
    down2up <- down2up+1
    }
    if (d[i-1]< 0 && d[i]<0){
    down2down <- down2down+1
    }

}

u2u <- up2up/(up2up+up2down)
u2d <- up2down/(up2up+up2down)
d2u <- down2up/(down2up+down2down)
d2d <- down2down/(down2up+down2down)

mattrans <- matrix(c(u2u, d2u, u2d, d2d),nrow=2,ncol=2)

colnames(mattrans) <- c("U","D")
row.names(mattrans) <- c("U","D")
rm(temp)
rm(sim)
temp<-0
sim <-data.frame(vector())

#now set states
state<- c(1,2,3,4,5)
nstate<-c(8000,9000,10000,11000,12000)

#current state
temp[1]<- 5

for (i in 1:100){
ins<- simula(mattrans,10)
 for (j in 2:20){
if (ins[j] == "U"){
    temp[j] <- temp[j-1]+1
}
else { temp[j] <- temp[j-1]-1}
}

sim <-rbind(sim, temp)
}

p1<-sum(rowSums(sim <=1))/(10*100)

p2<-sum(rowSums(sim == 2))/(10*100)

p3<-sum(rowSums(sim == 3))/(10*100)
p4<-sum(rowSums(sim == 4))/(10*100)
p5<-sum(rowSums(sim >= 5))/(10*100)

p <- c(p1,p2,p3,p4,p5)

cat("tansition matex is:\n")
matrix(mattrans, nrow = 2,ncol = 2)

f <- data.frame(nstate, p)
colnames(f)<-c("주가", "확률")

f[order(-f[,1]),]

cat("current state is:",temp[1], sep="\n")

2014년 4월 11일 금요일

markov chain for stock price

## NAME:  markovchain.py
##
## USAGE:  From shell prompt: python -i markovchain.py
##         or within interactive python environment, scrape and paste
## REQUIRED ARGUMENTS: none
##
## OPTIONS: none
## DESCRIPTION: Simulation of the cash management model as a Markov chain
## DIAGNOSTICS: None
## CONFIGURATION AND ENVIRONMENT: None
## DEPENDENCIES: None
## INCOMPATIBILITIES: None known
## PROVENANCE: Tue Aug 28, 2012  5:38 AM
## BUGS AND LIMITATIONS:  None Known
## FEATURES AND POTENTIAL IMPROVEMENTS:
## For the innermost for-loop, you can simply use
## newState = scipy.random.choice(1, 1, prob=T[state,])
## That shows more clearly what's going on
## On the other hand, if other languages don't have a function comparable to sample, then
## the current formulation is easier to adapt and make the algorithms similar.

## AUTHOR:  Steve Dunbar
## VERSION: 1.0 Tue Aug 28, 2012  5:38 AM
##          1.1 added some ouput formatting Mon Nov 12, 2012  6:03 AM
## KEYWORDS:   Markov chain, simulation, reflecting boundary
#! /usr/bin/env python

import scipy

n = 6;   # Top boundary, number of states 0..n is n+1
s = 4;   # Start and Reset state number 1 <= s <= n-1
p = 0.8;
steps = 1000;

T = scipy.diag( p*scipy.ones( n ), 1) + scipy.diag( (1-p)*scipy.ones( n ), -1);
T[0,1] = 0; T[0,s] = 1;
T[n,n-1] = 0; T[n,s] = 1;

# vector to hold the count of visits to each state during a cycle
count = scipy.zeros(n+1, float);
# Initialize the cyclelength
numberCycles = 0;
# Start in the state s
state = s;
# Make steps through the markov chain
for i in range(1, steps+1):
    x = 0;
    u = scipy.random.random(1);
    newState = state;
    for j in range(n+1):
    x = x + T[state, j];
    if (x >= u):
        newState = j;
        break;
    ## newState = scipy.random.choice(1, 1, prob=T[state,])
    state = newState;
    count[state] = count[state] + 1;
   
    if (state == n or state == 0):
        numberCycles = numberCycles + 1;
        avgVisits = count/numberCycles;
        avgCycleLength = i/numberCycles;

Wsk = avgVisits

print Wsk

theoreticalWsk = 2*(  s*(1-scipy.arange(0,n+1,dtype=float)/n) - scipy.maximum(s-scipy.arange(0,n+1), scipy.zeros(n+1, int) ) )

prob = [0] * (n+1)    # initialize probility as list data type
for i in range(1, n+1):
    prob[i] = Wsk[i]/sum(Wsk);

price = [2000, 2900, 3500, 4500, 5100, 5300, 5500]   # this is list not array
print "each price state in a cycle:\n", price
print "Average number of visits to each state in a cycle:\n ",  [round(n, 2) for n in prob]