Email         : ardifo.okta@student.matanauniversity.ac.id
RPubs       : https://rpubs.com/ardifo/
Jurusan : Statistika
Instagram : Ardifosyaa
Address : ARA Center, Matana University Tower
   Jl. CBD Barat Kav, RT.1, Curug Sangereng, Kelapa Dua, Tangerang, Banten 15810.

1 Question 1

Buatlah contoh kasus Optimasi Program Linear dalam Industri atau Bisnis!

Jawaban No. 1

Sebuah perusahaan PT. Maju Bersama memiliki sebuah pabrik yang akan memproduksi 2 jenis produk, yaitu Jaket Hoddie dan celana Jeans. Untuk memproduksi kedua produk diperlukan bahan baku benang, bahan baku kain jeans dan tenaga kerja. Maksimum penyediaan benang adalah 60 Kg perhari, kain jeans 30 kg per hari dan tenaga kerja 40 jam per bulan. Kebutuhan setiap unit produk akan bahan baku dan jam tenaga kerja dapat dilihat dalam tabel berikut ini:

library(readxl)
data_lopo <- read_excel("~/.A. KULIAH/Semester 7/Optimasi/Book1.xlsx")

## New names:
## • `` -> `...1`
## • `` -> `...3`
## • `` -> `...4`

print(data_lopo)

## # A tibble: 4 × 4
##   ...1                              `Kg bahan baku dan Jam tenaga …` ...3  ...4 
##   <chr>                             <chr>                            <chr> <chr>
## 1 Jenis bahan baku dan tenaga kerja Benang                           Kain… Maks…
## 2 Benang                            2                                3     60 kg
## 3 Kain Jeans                        0                                2     30 kg
## 4 Tenaga Kerja                      2                                1     40 j…

Kedua jenis produk memberikan keuntungan sebesar Rp40.000.000 untuk jenis benang dan Rp30.000.000 untuk jenis ki=ain jeans. Masalahnya yang kita punya dalah bagaimana menentukan jumlah unit setiap jenis produk yang akan diproduksi setiap hari agar keuntungan yagn diperoleh bisa maksimal.

Penyelesaian dengan Maksimasi Fungsi Linear Progamming

Menentukan Variabel terlebih dahulu

\(X_1 = Jenis Benang\)

\(X_2 = Kain Jeans\)

Fungsi Tujuan

\(Zmax = 40X_1 + 30X_2\)

library(lpSolve)
# Memasukan koefisien dari fungsi tujuan pada f.obj
f.obj <- c(40, 30)

Fungsi Kendala atau Batasan

\(2X_1 + 3X_2 <= 60 (benang)\)

\(2X_2 <= 30 (Kain Jeans)\)

\(2X_1 + X_2 <= 40 (tenaga kerja)\)

# memasukan koefisien fungsi kendala dalam bentuk matrix 
# dengan nrow menunjukan banyaknya kendala yaitu 3 dan angka yang diinput disusun perbaris sehingga by row = TRUE 
f.con <- matrix(c(2, 3,
                  0, 2,
                  2, 1), nrow = 3, byrow= TRUE)

# Memasukkan tanda pertidaksamaan pada setiap kendala
f.dir <- c("<=",
           "<=",
           "<=")

# Memasukkan koefisien ruas kanan
f.rhs <- c(60,
           30,
           40)

Mendapatkan solusi Optimal

# Keuntungan Maksimum

lp("max", f.obj, f.con, f.dir, f.rhs)

## Success: the objective function is 900

# Nilai Variabel agar mencapai keuntungan maksimum
lp("max", f.obj, f.con, f.dir, f.rhs)$solution

## [1] 15 10

# Koefisien sensitivitas
lp("max", f.obj, f.con, f.dir, f.rhs, compute.sens=TRUE)$sens.coef.from

## [1] 20 20

lp("max", f.obj, f.con, f.dir, f.rhs, compute.sens=TRUE)$sens.coef.to

## [1] 60 60

# Bayangan/harga ganda dari kendala dan variabel keputusan
lp("max", f.obj, f.con, f.dir, f.rhs, compute.sens=TRUE)$duals

## [1]  5  0 15  0  0

# Batas bawah batasan harga bayangan/ganda dan masing masing variabel keputusan
lp("max", f.obj, f.con, f.dir, f.rhs, compute.sens=TRUE)$duals.from

## [1]  4e+01 -1e+30  3e+01 -1e+30 -1e+30

# Batas atas batasan harga bayangan/ganda dan masing masing variabel keputusan
lp("max", f.obj, f.con, f.dir, f.rhs, compute.sens=TRUE)$duals.to

## [1] 7e+01 1e+30 6e+01 1e+30 1e+30

Kesimpulan

Nilai z maksimum yang dapat diperoleh saat memenuhi kendala yang diberikan adalah 900

(keuntungan sebesar Rp 900 juta)

di mana

\(X_2=15\)

\(X_2=10\)

Koefisien sensitivitas berubah dari 20 dan 20 menjadi 60 dan 60

Bayangan/harga ganda dari kendala adalah 5 , 0 , 15

variabel keputusan masing-masing adalah 0 dan 0

Batas bawah harga bayangan adalah 4e+01, -1e+30, dan 3e+01

variabel keputusan berturut-turut adalah -1e+30 dan -1e+30

Batas atas harga bayangan adalah 17e+01, 1e+30, dan 6e+01

variabel keputusan masing-masing adalah 1e+30 dan 1e+30

2 Question2

Buatlah contoh kasus Optimasi Program NonLinear dalam Industri atau Bisnis!

Ketika hubungan antar variabel tidak linier, Anda dapat mencoba:

mentransformasikan data untuk melinearisasi hubungan
- paskan fungsi non-linier ke data (gunakan contoh nls)
- paskan model polinomial atau spline ke data (gunakan contoh paket growthrates)

Regresi Linear : variabel terikat = konstanta + parameter x variabel independen + p x IV +…. \(y=B_0+B_1X_1+B_2X_2+…\)

Ini tidak berarti Anda tidak bisa menyesuaikan kurva! Yang membuatnya linier adalah parameternya yang linier. Anda mungkin memiliki variabel independen nonlinier misalnya) \(y=B_0+B_1X_1+B_2X_2\). Daripada mengubah variabel independen secara eksponensial, Anda dapat menggunakan suku log atau invers, dll.

regresi nonlinier : Ada lagi. Modelnya bisa seperti: \(y=B_1xcos(X+B_4)+B_2xcos(2∗X+B_4)+B_3\). Ini membuatnya penting bagi Anda untuk melakukan penelitian untuk memahami bentuk fungsional apa yang mungkin data anda ambil.

Pendekatan Non Linear Least Squares

Kuadrat terkecil nonlinier adalah cara yang baik untuk memperkirakan parameter agar sesuai dengan data nonlinier. Menggunakan fungsi linier untuk memperkirakan fungsi nonlinier dan secara berulang bekerja untuk menemukan nilai parameter terbaik.

# simulasi beberapa data
# set generator nomor acak

set.seed(2019)
x<-seq(0,50,1)
y<-((runif(1,10,20)*x)/(runif(1,0,10)+x))+rnorm(51,0,1)

# model sederhana nls untuk  menemukan nilai awal yang baik bagi parameter meskipun itu melontarkan peringatan

# Persamaan Michaelis-Menten
m <-nls (y ~ a * x / (b + x))

## Warning in nls(y ~ a * x/(b + x)): No starting values specified for some parameters.
## Initializing 'a', 'b' to '1.'.
## Consider specifying 'start' or using a selfStart model

## Nonlinear regression model
##   model: y ~ a * x/(b + x)
##    data: parent.frame()
##      a      b 
## 17.908  7.645 
##  residual sum-of-squares: 51.78
## 
## Number of iterations to convergence: 6 
## Achieved convergence tolerance: 1.077e-06

# Dapatkan beberapa perkiraan tentang kesesuaian
cor(y,predict(m))

## [1] 0.9658715

# Plot
plot(x,y)
lines(x,predict(m),col="blue",lwd=3)

# Mensimulasikan beberapa data, ini tanpa pengetahuan apriori tentang nilai parameter
x <- seq(0,50,1)
y <- runif(1,5,15)*exp(-runif(1,0.01,0.05)*x)+rnorm(51,0,0.5)

# Perkirakan secara visual beberapa nilai parameter awal
plot(x,y)

# Dari grafik ini mengatur perkiraan nilai awal
a_start<-8               # param a adalah nilai y saat x = 0
b_start<-log(0.1)/(50*8) # b adalah tingkat kerusakan. k = log (A) / (A (initial) * t)

# Model
m1<-nls(y~a*exp(-b*x),start=list(a=a_start,b=b_start))
m1

## Nonlinear regression model
##   model: y ~ a * exp(-b * x)
##    data: parent.frame()
##       a       b 
## 8.59328 0.03544 
##  residual sum-of-squares: 7.588
## 
## Number of iterations to convergence: 6 
## Achieved convergence tolerance: 7.615e-07

# Mendapatkan beberapa perkiraan tentang kesesuaian (harus mendekati 1)
cor(y,predict(m1))

## [1] 0.9833285

# Plot the fit
plot(x,y)
lines(x,predict(m1),col="blue",lty=2,lwd=3)

library(deSolve) # Paket ini bagus untuk menyelesaikan persamaan diferensial

# Mensimulasikan beberapa pertumbuhan populasi dari persamaan logistik dan mengestimasi parameter menggunakan nls

# Mendefinisikan suatu fungsi
log_growth <- function(Time, State, Pars) {
  with(as.list(c(State, Pars)), {
    dN <- R*N*(1-N/K)
    return(list(c(dN)))
  })
}

# Waktu adalah interval waktu, Status adalah nama variabel, Pars, dan parameter

# Pars, N_int, times digunakan untuk mensimulasikan data
pars  <- c(R=0.2,K=1000) # Parameter untuk logisitc growth
N_ini <- c(N=1) # Inisial Nomor
times <- seq(0, 50, by = 1) # Langkah waktu untuk mengevaluasi

# ODE (ordinary differential equation)
out <- ode(N_ini, times, log_growth, pars)
plot(out)

# Tambahkan beberapa variasi acak padanya
N_obs<-out[,2]+rnorm(51,0,50)
# Hapus angka kurang dari 1
N_obs<-ifelse(N_obs<1,1,N_obs)
# plot
plot(times,N_obs)

# Tidak memiliki nilai awal hanya berfungsi terkadang dengan data dan fungsi sederhana seperti pada contoh pertama.

# Perhatikan bagaimana m3 TIDAK AKAN bekerja. Ingat nls secara berulang berjalan untuk berkumpul dan ini tidak akan berkumpul di bawah 30 iterasi
#m3<-nls(N_obs~K*N0*exp(R*times)/(K+N0*(exp(R*times)-1)))

# Mendapatkan Initial memberikan tebakan pada nilai parameter berdasarkan data

# SSlogis adalah model yang memulai sendiri
SS<-getInitial(N_obs~SSlogis(times,alpha,xmid,scale),data=data.frame(N_obs=N_obs,times=times))
# Ikuti persamaan berikut: N(t)=alpha/(1+e+((xmid-t)/scale)
SS

##      alpha       xmid      scale 
## 991.798725  34.112932   5.045119

# Perlu melakukan beberapa aljabar untuk mendapatkan parameterisasi yang benar
# Kita menggunakan parameterisasi yang berbeda
K_start<-SS["alpha"]
R_start<-1/SS["scale"]
N0_start<-SS["alpha"]/(exp(SS["xmid"]/SS["scale"])+1)

# Rumus (tidak disiapkan sebagai persamaan diferensial) untuk model
log_formula<-formula(N_obs~K*N0*exp(R*times)/(K+N0*(exp(R*times)-1)))

# Cocokan dengan modelnya
m4<-nls(log_formula,start=list(K=K_start,R=R_start,N0=N0_start))

# Estimasi parameter
summary(m4)

## 
## Formula: N_obs ~ K * N0 * exp(R * times)/(K + N0 * (exp(R * times) - 1))
## 
## Parameters:
##           Estimate Std. Error t value Pr(>|t|)    
## K.alpha  991.79872   29.86851  33.205   <2e-16 ***
## R.scale    0.19821    0.01338  14.817   <2e-16 ***
## N0.alpha   1.14659    0.47931   2.392   0.0207 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 44.08 on 48 degrees of freedom
## 
## Number of iterations to convergence: 0 
## Achieved convergence tolerance: 4.177e-06

# Dapatkan beberapa estimasi kebaikan kecocokan (harus lebih dekat dengan 1)
cor(N_obs,predict(m4))

## [1] 0.9925536

# Plot
plot(times,N_obs)
lines(times,predict(m4),col="darkred",lty=2,lwd=3)

Pendekatan kemungkinan maksimum lihat paket nlme. Metode yang sedikit lebih kuat dan andal daripada nls.

Growthrates package

Menemukan yang paling cocok agak mengganggu karena menyesuaikan fungsi yang agak sederhana. Pasang spline! Spline tidak hanya membuat garis melalui data, tetapi juga melewati setiap titik data dan menyesuaikan ploynomial kubik antara dua titik. Jadi spline adalah snapshot yang paling cocok. Memungkinkan interpolasi yang baik. fungsi fit_spline menyesuaikan spline ke satu dan menyelesaikan sistem persamaan kubik.

library(growthrates)

## Loading required package: lattice

time<-rep(1:25, 8)
y<-grow_logistic(time, c(y0=0.2, mumax=0.3, K=4))[,"y"]
y<-jitter(y, factor=1, amount=0.3)
data<-data.frame(cbind(time,y))
data<-data[order(data$y),]
data$ID<-c(rep(1:10,20))
plot(time,y)

splitted.data <- multisplit(data, c("ID"))

# Menunjukkan eksperimen mana yang ada di splitted.data
names(splitted.data)

##  [1] "1"  "2"  "3"  "4"  "5"  "6"  "7"  "8"  "9"  "10"

# Dapatkan tabel dari satu individu, atau bisa menjadi percobaan, atau kelompok, atau blok, dll.
dat <- splitted.data[["1"]]
dat[1, 2] = 0.03

fit0 <- fit_spline(dat$time, dat$y)
plot(fit0, col="red", lwd=3)

summary(fit0)

## Fitted smoothing spline:
## Call:
## smooth.spline(x = time, y = ylog)
## 
## Smoothing Parameter  spar= 0.6721364  lambda= 0.006882113 (12 iterations)
## Equivalent Degrees of Freedom (Df): 3.685259
## Penalized Criterion (RSS): 2.029662
## GCV: 0.3015
## 
## Parameter values of exponential growth curve:
## Maximum growth at x= 1.001382 , y= 0.1644875 
## y0 = 0.1145593 
## mumax = 0.3612431 
## 
## r2 of log transformed data= 0.8652797

# Bisa mendapatkan tingkat pertumbuhan maksimum dari sini. Plot akan menampilkan kurva eksponensial (garis biru) dari nilai maks
# Inisial Parameter
plot(fit0, col="blue", lwd=3)
p <- c(coef(fit0), K = max(dat$y))

# Hindari parameter negatif
lower = c(y0 = 0, mumax = 0, K = 0)

fit3<-fit_growthmodel(grow_logistic, p=p, time=data$time, y=data$y)

## Warning in regularize.values(x, y, ties, missing(ties), na.rm = na.rm):
## collapsing to unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties), na.rm = na.rm):
## collapsing to unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties), na.rm = na.rm):
## collapsing to unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties), na.rm = na.rm):
## collapsing to unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties), na.rm = na.rm):
## collapsing to unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties), na.rm = na.rm):
## collapsing to unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties), na.rm = na.rm):
## collapsing to unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties), na.rm = na.rm):
## collapsing to unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties), na.rm = na.rm):
## collapsing to unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties), na.rm = na.rm):
## collapsing to unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties), na.rm = na.rm):
## collapsing to unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties), na.rm = na.rm):
## collapsing to unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties), na.rm = na.rm):
## collapsing to unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties), na.rm = na.rm):
## collapsing to unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties), na.rm = na.rm):
## collapsing to unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties), na.rm = na.rm):
## collapsing to unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties), na.rm = na.rm):
## collapsing to unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties), na.rm = na.rm):
## collapsing to unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties), na.rm = na.rm):
## collapsing to unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties), na.rm = na.rm):
## collapsing to unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties), na.rm = na.rm):
## collapsing to unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties), na.rm = na.rm):
## collapsing to unique 'x' values

lines(fit3, col="red", lwd=3)

summary(fit3)

## 
## Parameters:
##       Estimate Std. Error t value Pr(>|t|)    
## y0    0.188670   0.012562   15.02   <2e-16 ***
## mumax 0.301080   0.007458   40.37   <2e-16 ***
## K     4.027982   0.029659  135.81   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1713 on 197 degrees of freedom
## 
## Parameter correlation:
##            y0   mumax       K
## y0     1.0000 -0.9448  0.5116
## mumax -0.9448  1.0000 -0.6615
## K      0.5116 -0.6615  1.0000

# Dapatkan parameter untuk seluruh populasi

fit4<-all_growthmodels(
          y ~ grow_logistic(time, parms) | ID,
          data = data, p = p, lower = lower) 
# Memberikan parameter untuk setiap individu

results(fit4, extended=TRUE) # Diperpanjang memberi Anda perkiraan waktu nilai saturasi.

##    ID        y0     mumax        K turnpoint     sat1     sat2     sat3
## 1   1 0.1775788 0.3186817 3.922830  9.567004 13.69952 15.84286 18.80644
## 2   2 0.2131942 0.2973178 3.962413  9.643233 14.07267 16.37004 19.54655
## 3   3 0.1958260 0.2873056 4.118416 10.432393 15.01620 17.39362 20.68083
## 4   4 0.1639432 0.3095515 4.076656 10.248580 14.50299 16.70954 19.76054
## 5   5 0.1514473 0.3280776 3.963867  9.832374 13.84654 15.92849 18.80720
## 6   6 0.1988192 0.2916016 4.077459 10.188041 14.70430 17.04671 20.28550
## 7   7 0.2349602 0.2860918 4.040224  9.733663 14.33695 16.72443 20.02562
## 8   8 0.1929363 0.2827270 4.130590 10.667467 15.32552 17.74143 21.08189
## 9   9 0.1830012 0.3046824 4.035454 10.000489 14.32289 16.56470 19.66445
## 10 10 0.1680191 0.3139180 4.027611  9.984249 14.17949 16.35534 19.36390
##           r2
## 1  0.9896373
## 2  0.9803218
## 3  0.9865337
## 4  0.9893162
## 5  0.9840925
## 6  0.9897703
## 7  0.9846571
## 8  0.9887043
## 9  0.9836142
## 10 0.9825879

Bagaimana ID diatur sekarang, Anda akan mendapatkan estimasi parameter untuk setiap individu. Dapat mengubah ID menjadi perlakuan, percobaan, pemblokiran, dll. Berdasarkan data Anda. Jika Anda menggantinya dengan perawatan, Anda akan mendapatkan parameter untuk setiap perawatan. Dapat menyertakan lebih banyak dengan menambahkan tanda tambah di antara mereka.

Parameter: y0 = nilai y awal; mumax = tingkat pertumbuhan intrinsik; K = ukuran asimtotik

Waktu saturasi: titik balik = waktu titik balik (saturasi 50%); sat1 = waktu minimum turunan ke-2 (laju pertumbuhan minimum); sat2 = waktu intersep antara kenaikan paling tajam (garis singgung mumaks) dan K; sat3 = waktu ketika ukuran asimtotik tercapai

Ini menunjukkan bagaimana menyesuaikan persamaan logistik, tetapi paket ini cocok dengan banyak kurva yang berbeda! Jika Anda lebih tertarik untuk mendefinisikan forumlas yang lebih kompleks dan mendefinisikan fungsi log-liklihood mereka, saya merekomendasikan paket bblme dan menggunakan fungsi mle2

3 Question 3

Buatlah contoh kasus Optimasi Portofolio dengan memilih 5 saham terbaik di Indonesia dengan cara melakukan analisis keuangan terlebih dahulu pada saham LQ45!

Dalam sebuah index Indonesia mempunyai 45 perusahaan yang sudah IPO, perusahaan-perusahaan tersebut telah membuka saham di pasar saham melalaui Bursa Efek Indonesia. Seorang anak muda berinvestasi saham pada perusahaan -perusahaan yag sudah masuk di daftar LQ45 tetapi ia bingung sebab dari sekian banyak perusahaan tersebut ia harus memilih 5 perusahaan untuk berinvestasi. Agar investasi nya menghasilkan profit dan menghindarkan loss maka kita akan melakukan optimsi dengan memilih 5 saham terbaik untuk tempat berinvestasi dengan menganalisakeuangan terlebih daulu. berikut ini adalah kelima saham yang di rekomendasikan untuk berinvestasi yetaapi sebelum berinvestasi maka dilakukan analisis terlebih dahulu.

Data Saham Data saham yang saya gunakan ialah sebagai berikut.

Kode | Nama Stok | Sektor

TPIA | Chandra Asri Petrochemical Tbk. |Basic Industry and Chemicals
ANTM |Aneka Tambang Tbk. [S] |Metal and Mining, 23
LPKR | Lippo Karawaci Tbk. [S] |Property and Real Estate, 61
KLBF | Kalbe Farma Tbk. [S] |Pharmaceuticals, 300
INKP | Indah Kiat Pulp & Paper |Basic Industry and Chemicals

library(tidyquant)

## Loading required package: lubridate

## 
## Attaching package: 'lubridate'

## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union

## Loading required package: PerformanceAnalytics

## Loading required package: xts

## Loading required package: zoo

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

## 
## Attaching package: 'PerformanceAnalytics'

## The following object is masked from 'package:graphics':
## 
##     legend

## Loading required package: quantmod

## Loading required package: TTR

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

library(plotly)

## Loading required package: ggplot2

## 
## Attaching package: 'plotly'

## The following object is masked from 'package:ggplot2':
## 
##     last_plot

## The following object is masked from 'package:stats':
## 
##     filter

## The following object is masked from 'package:graphics':
## 
##     layout

library(timetk)

#daftar data saham yang ingin diimport
tick <- c('TPIA.JK', 'ANTM.JK', 'LPKR.JK', 'KLBF.JK', 'INKP.JK')

price_data <- tq_get(tick,
                     from = '2020-01-01',
                     to = Sys.Date(),
                     get = 'stock.prices')
head(price_data)

## # A tibble: 6 × 8
##   symbol  date        open  high   low close   volume adjusted
##   <chr>   <date>     <dbl> <dbl> <dbl> <dbl>    <dbl>    <dbl>
## 1 TPIA.JK 2020-01-02 2612. 2619. 2544. 2569. 18964400    2545.
## 2 TPIA.JK 2020-01-03 2569. 2569. 2419. 2569. 18256800    2545.
## 3 TPIA.JK 2020-01-06 2550  2550  2431. 2438. 40981200    2415.
## 4 TPIA.JK 2020-01-07 2462. 2488. 2375  2406. 19280400    2384.
## 5 TPIA.JK 2020-01-08 2394. 2400  2319. 2325  24982800    2304.
## 6 TPIA.JK 2020-01-09 2325  2381. 2325  2369. 20129200    2347.

Pengembalian Menghitung return harian untuk saham-saham diatas dengan menggunakan pengembalian logaritmik.

log_ret_tidy <- price_data %>%
  group_by(symbol) %>%
  tq_transmute(select = adjusted,
               mutate_fun = periodReturn,
               period = 'daily',
               col_rename = 'ret',
               type = 'log')

#log_ret_tidy$ret <-round(log_ret_tidy$ret,4)

library(DT)
datatable(log_ret_tidy)

selanjutnya kita menggunakan fungsispread() untu mengubah data tabel menjadi format yang agak lebar menggunakan xts() untuk mengubahnya menjadi objek time series.

library(tidyr)

log_ret_xts <- log_ret_tidy %>%
  spread(symbol, value = ret) %>%
  tk_xts()

## Warning: Non-numeric columns being dropped: date

## Using column `date` for date_var.

datatable(log_ret_xts)

Rata-rata Pengembalia selanjutnya menghitung rata-rata pengembalian harian untuk setiap aset dari beberapa perusahaan tersebut.

mean_ret <- colMeans(log_ret_xts)
print(round(mean_ret,4))

## ANTM.JK INKP.JK KLBF.JK LPKR.JK TPIA.JK 
##  0.0013  0.0003  0.0004 -0.0015  0.0000

Matriks Kovariansi Untuk menghitung matriks kovarians untuk semua stok dengan mengalikannya dengan 360 dalam format tahunan.

cov_mat <- cov(log_ret_xts) * 360
print(round(cov_mat,4))

##         ANTM.JK INKP.JK KLBF.JK LPKR.JK TPIA.JK
## ANTM.JK  0.4809  0.1922  0.0603  0.1121  0.0450
## INKP.JK  0.1922  0.3939  0.0689  0.0962  0.0575
## KLBF.JK  0.0603  0.0689  0.1900  0.0562  0.0216
## LPKR.JK  0.1121  0.0962  0.0562  0.3456  0.0360
## TPIA.JK  0.0450  0.0575  0.0216  0.0360  0.2058

Penerapan Metode Portofolio

Langkah-langkah untuk memperhitungkan pengembalian dan risiko portofolio sebagai berikut:

Bobok acak
rata-rata pengembalian aset
Risiko Portofolio
Bobot Portofolio

wts <- runif(n = length(tick))
wts <- wts/sum(wts)

port_returns <- (sum(wts * mean_ret) + 1)^252 - 1
port_risk <- sqrt(t(wts) %*% (cov_mat %*% wts))
sharpe_ratio <- port_returns/port_risk

Selanjutnya melakukan proses pembentukan portofolio secara acak dengan simuasi 5000 kali untuk memastikan signifikansinya secara statistik.

num_port <- 5000

# matriks kosong untuk menyimpan bobot
all_wts <- matrix(nrow = num_port, ncol = length(tick))

# vektor kosong untuk menyimpan pengembalian portofolio
port_returns <- vector('numeric', length = num_port)

# vektor kosong untuk menyimpan standar deviasi portofolio
port_risk <- vector('numeric', length = num_port)

#vektor kosong untuk menyimpan portofolio sharpe ratio
sharpe_ratio <- vector('numeric', length = num_port)

Selanjutnya akan dilakukan loop 5000 kali.

for(i in seq_along(port_returns)) {
  
  wts<- runif(length(tick))
  wts <- wts/sum(wts)
  
  #menyimpan proporsi saham dalam matriks
  all_wts[i,] <- wts
  
  #pengembalian portofolio
  port_ret <- sum(wts * mean_ret)
  port_ret <- ((port_ret +1)^360) - 1
  
  #menyimpan nilai pengembalian portofolio
  port_returns[i] <- port_ret
  
  # membuat dan menyimpan risiko portofolio
  port_sd <- sqrt(t(wts) %*% (cov_mat %*% wts))
  port_risk[i] <-port_sd
  
  #Membuat dan menyimpan portofolio sharpe ratio
  # dengan asumsi tingkat bebas risiko 0%
  sr <- port_ret/port_sd
  sharpe_ratio[i] <- sr
  
}

Selanjutnya, membuat tabel data untuk menyimpan semua nilai secara bersamaan.

# menyimpan nilai dalam tabel
portfolio_values <- tibble(Return  = port_returns,
                            Risk = port_risk,
                            SharpeRatio = sharpe_ratio)

# mengubah matriks mwnjadi tibble
all_wts <- tk_tbl(all_wts)

## Warning in tk_tbl.data.frame(as.data.frame(data), preserve_index,
## rename_index, : Warning: No index to preserve. Object otherwise converted to
## tibble successfully.

# mengubah nama kolom
colnames(all_wts) <- colnames(log_ret_xts)

# menyisir (susun) semua nilai bersama-sama
portfolio_values <- tk_tbl(cbind(all_wts, portfolio_values))

## Warning in tk_tbl.data.frame(cbind(all_wts, portfolio_values)): Warning: No
## index to preserve. Object otherwise converted to tibble successfully.

#print hasil dalam tabel
datatable(portfolio_values)

Karena sudah memiliki bobot di setiap aset dengan risiko dan pengembalian bersama dengan sharpe ratio dari setiap portofolio mari lihat portofolio yang paling penting dengan varians minimum dan juga sharpe ratio tertinggi.

Varians Minimum

library(forcats)

min_var <- portfolio_values[which.min(portfolio_values$Risk),]

p <- min_var %>%
  gather(ANTM.JK:KLBF.JK, key = Asset,
         value = Weights) %>%
  mutate(Asset = as.factor(Asset)) %>%
  ggplot(aes(x = fct_reorder(Asset,Weights), y = Weights, fill = Asset)) +
  geom_bar(stat = 'identity') +
  theme_minimal() +
  labs(x = 'Aset', 
       y = 'Bobot', 
       title = "Bobot Portofolio dengan Variansi Minimum") +
  scale_y_continuous(labels = scales::percent) +
  theme(legend.position = "none" )

ggplotly(p)

Setelah kita amati berdasarkan dengan hasil pemograman di atas, kita mendapati mayoritas portofolio diinvestasikan pada saham Kalbe FArma dan Aneka Tambang.

Portofolio Tangensi

# maksimum sharpe ratio
max_sr <- portfolio_values[which.max(portfolio_values$SharpeRatio),]

rii <- max_sr %>%
  gather(ANTM.JK:KLBF.JK, key = Asset,
         value = Weights) %>%
  mutate(Asset = as.factor(Asset)) %>%
  ggplot(aes(x = fct_reorder(Asset, Weights), y = Weights, fill = Asset)) +
  geom_bar(stat = 'identity') +
  theme_minimal() +
  labs( x = 'Aset',
        y = 'Bobot',
        title = "Bobot Portfolio Tangensi(Max Sharpe Ratio)") +
  scale_y_continuous(labels = scales::percent) +
  theme(legend.position = "none")

ggplotly(rii)

Portofolio dengan Ratio Sharpe tertinggi hanya pada perusahaan Aneka Tambang dan Kalbe Farma. Batas Efisien Portofolio

p <- portfolio_values %>%
  ggplot(aes(x = Risk, y = Return, color = SharpeRatio)) +
  geom_point() +
  theme_classic() +
  scale_y_continuous(labels = scales::percent) +
  scale_x_continuous(labels = scales::percent) +
  labs(x = 'Risiko Tahunan',
       y = 'Pengembalian Tahunan',
       title = "Optimasi Portofolio & Perbatasan yang Efisien") +
  geom_point(aes(x = Risk,y = Return), data = min_var, color = 'red') +
  geom_point(aes(x = Risk,y = Return), data = max_sr, color = 'red') +
  annotate('text', x = 0.71, y = 0.81, label = "Portofolio Tangensi") +
  annotate('text', x = 0.290, y = 0.250, label = "Portofolio Varians minimum") +
  annotate(geom = 'segment', x = 0.6920, xend = 0.6920,  y = 0.500, 
           yend = 0.6800, color = 'red', arrow = arrow(type = "open")) +
  annotate(geom = 'segment', x = 0.2870, xend = 0.2870,  y = 0.06, 
           yend = 0.200, color = 'red', arrow = arrow(type = "open"))
  

ggplotly(p)

FINAL EXAM

Optimasi

Ardifo Okta Yudistio (20194920001)

December 22, 2022

1 Question 1

2 Question2

3 Question 3