#input library
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(ggplot2)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(data.table)
## 
## Attaching package: 'data.table'
## The following objects are masked from 'package:dplyr':
## 
##     between, first, last

1 Input Data dan visualisasi Data

# Membuat data frame di R
data <- data.frame(
  Tahun = c(1980:2023),
   Produktivitas = c(4169.442, 3951.818, 3502.381, 3750.342, 3528.127, 3344.306, 3815.466, 4041.731,
                    3812.039, 3869.066, 3858.184, 3824.65, 3853.345, 3824.334, 3947.745, 3921.171,
                    3961.818, 3661.605, 4460.42, 4716.56, 4398.982, 4334.105, 4970.317, 5195.834,
                    4964.834, 5101.404, 5212.358, 5220.077, 5389.627, 5391.495, 5675.321, 5176.68,
                    5595.152, 5443.138, 5232.585, 5198.858, 5325.21, 5796.57, 5819.193, 5863.358,
                    2305.98, 3483.08, 6012.53, 5967.92),
  Suhu = c(26, 26.6, 26.4, 26.7, 26.1, 26.3, 26.3, 26.6, 26.5, 26.4, 26.6, 26.4, 26.5, 26.5, 26.4,
           26.6, 26.5, 26.6, 27.1, 26.4, 26.5, 26.7, 26.9, 26.8, 26.8, 26.9, 26.8, 26.8, 26.6, 27.0,
           27.1, 26.8, 26.9, 27.0, 27.0, 27.1, 27.4, 27.0, 27.1, 27.2, 27.0, 27.3, 27.0, 27.1),
  CurahHujan = c(80.4, 78.3, 77.8, 78.6, 75.8, 78.4, 81.0, 79.7, 81.1, 79.3, 79.3, 79.8, 78.2, 76.6,
                 81.5, 78.2, 75.4, 76.4, 74.2, 80.4, 81.2, 82.5, 80.2, 75.9, 79.2, 81.0, 80.9, 81.2,
                 30.5, 80.3, 50.8, 60.2, 83.2, 50.5, 70.6, 70.9, 90.4, 80.6, 81.3, 90.9, 72.2, 70.6,
                 30.5, 80.3)
)

# Tampilkan data frame
print(data)
##    Tahun Produktivitas Suhu CurahHujan
## 1   1980      4169.442 26.0       80.4
## 2   1981      3951.818 26.6       78.3
## 3   1982      3502.381 26.4       77.8
## 4   1983      3750.342 26.7       78.6
## 5   1984      3528.127 26.1       75.8
## 6   1985      3344.306 26.3       78.4
## 7   1986      3815.466 26.3       81.0
## 8   1987      4041.731 26.6       79.7
## 9   1988      3812.039 26.5       81.1
## 10  1989      3869.066 26.4       79.3
## 11  1990      3858.184 26.6       79.3
## 12  1991      3824.650 26.4       79.8
## 13  1992      3853.345 26.5       78.2
## 14  1993      3824.334 26.5       76.6
## 15  1994      3947.745 26.4       81.5
## 16  1995      3921.171 26.6       78.2
## 17  1996      3961.818 26.5       75.4
## 18  1997      3661.605 26.6       76.4
## 19  1998      4460.420 27.1       74.2
## 20  1999      4716.560 26.4       80.4
## 21  2000      4398.982 26.5       81.2
## 22  2001      4334.105 26.7       82.5
## 23  2002      4970.317 26.9       80.2
## 24  2003      5195.834 26.8       75.9
## 25  2004      4964.834 26.8       79.2
## 26  2005      5101.404 26.9       81.0
## 27  2006      5212.358 26.8       80.9
## 28  2007      5220.077 26.8       81.2
## 29  2008      5389.627 26.6       30.5
## 30  2009      5391.495 27.0       80.3
## 31  2010      5675.321 27.1       50.8
## 32  2011      5176.680 26.8       60.2
## 33  2012      5595.152 26.9       83.2
## 34  2013      5443.138 27.0       50.5
## 35  2014      5232.585 27.0       70.6
## 36  2015      5198.858 27.1       70.9
## 37  2016      5325.210 27.4       90.4
## 38  2017      5796.570 27.0       80.6
## 39  2018      5819.193 27.1       81.3
## 40  2019      5863.358 27.2       90.9
## 41  2020      2305.980 27.0       72.2
## 42  2021      3483.080 27.3       70.6
## 43  2022      6012.530 27.0       30.5
## 44  2023      5967.920 27.1       80.3
ggplot(data, aes(x = Tahun, y = Produktivitas)) +
  geom_line(color = "blue") +
  ggtitle("Tren Produktivitas Pertanian (1980-2023)") +
  xlab("Tahun") +
  ylab("Produktivitas (ton/ha)") +
  theme_minimal()

ggplot(data, aes(x = Tahun, y = Suhu)) +
  geom_line(color = "blue") +
  ggtitle("Tren Suhu (1980-2023)") +
  xlab("Tahun") +
  ylab("Produktivitas (ton/ha)") +
  theme_minimal()

ggplot(data, aes(x = Tahun, y = CurahHujan)) +
  geom_line(color = "blue") +
  ggtitle("Tren Curah Hujan (1980-2023)") +
  xlab("Tahun") +
  ylab("Produktivitas (ton/ha)") +
  theme_minimal()

2 Membuat data menjadi logaritma natural (ln)

data$lnProduktivitas <- log(data$Produktivitas)
data$lnSuhu <- log(data$Suhu)
data$lnCurahHujan <- log(data$CurahHujan)

data
##    Tahun Produktivitas Suhu CurahHujan lnProduktivitas   lnSuhu lnCurahHujan
## 1   1980      4169.442 26.0       80.4        8.335537 3.258097     4.387014
## 2   1981      3951.818 26.6       78.3        8.281931 3.280911     4.360548
## 3   1982      3502.381 26.4       77.8        8.161198 3.273364     4.354141
## 4   1983      3750.342 26.7       78.6        8.229602 3.284664     4.364372
## 5   1984      3528.127 26.1       75.8        8.168522 3.261935     4.328098
## 6   1985      3344.306 26.3       78.4        8.115014 3.269569     4.361824
## 7   1986      3815.466 26.3       81.0        8.246818 3.269569     4.394449
## 8   1987      4041.731 26.6       79.7        8.304428 3.280911     4.378270
## 9   1988      3812.039 26.5       81.1        8.245919 3.277145     4.395683
## 10  1989      3869.066 26.4       79.3        8.260768 3.273364     4.373238
## 11  1990      3858.184 26.6       79.3        8.257952 3.280911     4.373238
## 12  1991      3824.650 26.4       79.8        8.249222 3.273364     4.379524
## 13  1992      3853.345 26.5       78.2        8.256697 3.277145     4.359270
## 14  1993      3824.334 26.5       76.6        8.249140 3.277145     4.338597
## 15  1994      3947.745 26.4       81.5        8.280900 3.273364     4.400603
## 16  1995      3921.171 26.6       78.2        8.274146 3.280911     4.359270
## 17  1996      3961.818 26.5       75.4        8.284458 3.277145     4.322807
## 18  1997      3661.605 26.6       76.4        8.205657 3.280911     4.335983
## 19  1998      4460.420 27.1       74.2        8.402998 3.299534     4.306764
## 20  1999      4716.560 26.4       80.4        8.458835 3.273364     4.387014
## 21  2000      4398.982 26.5       81.2        8.389128 3.277145     4.396915
## 22  2001      4334.105 26.7       82.5        8.374270 3.284664     4.412798
## 23  2002      4970.317 26.9       80.2        8.511239 3.292126     4.384524
## 24  2003      5195.834 26.8       75.9        8.555612 3.288402     4.329417
## 25  2004      4964.834 26.8       79.2        8.510135 3.288402     4.371976
## 26  2005      5101.404 26.9       81.0        8.537271 3.292126     4.394449
## 27  2006      5212.358 26.8       80.9        8.558788 3.288402     4.393214
## 28  2007      5220.077 26.8       81.2        8.560267 3.288402     4.396915
## 29  2008      5389.627 26.6       30.5        8.592231 3.280911     3.417727
## 30  2009      5391.495 27.0       80.3        8.592578 3.295837     4.385770
## 31  2010      5675.321 27.1       50.8        8.643882 3.299534     3.927896
## 32  2011      5176.680 26.8       60.2        8.551919 3.288402     4.097672
## 33  2012      5595.152 26.9       83.2        8.629656 3.292126     4.421247
## 34  2013      5443.138 27.0       50.5        8.602111 3.295837     3.921973
## 35  2014      5232.585 27.0       70.6        8.562661 3.295837     4.257030
## 36  2015      5198.858 27.1       70.9        8.556194 3.299534     4.261270
## 37  2016      5325.210 27.4       90.4        8.580207 3.310543     4.504244
## 38  2017      5796.570 27.0       80.6        8.665022 3.295837     4.389499
## 39  2018      5819.193 27.1       81.3        8.668917 3.299534     4.398146
## 40  2019      5863.358 27.2       90.9        8.676478 3.303217     4.509760
## 41  2020      2305.980 27.0       72.2        7.743261 3.295837     4.279440
## 42  2021      3483.080 27.3       70.6        8.155672 3.306887     4.257030
## 43  2022      6012.530 27.0       30.5        8.701601 3.295837     3.417727
## 44  2023      5967.920 27.1       80.3        8.694154 3.299534     4.385770

3 Model Ekspetasi Cobb-Douglas

model_ekspektasi <- lm(lnProduktivitas ~ lnSuhu + lnCurahHujan, data=data)

summary(model_ekspektasi)
## 
## Call:
## lm(formula = lnProduktivitas ~ lnSuhu + lnCurahHujan, data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.75263 -0.06896  0.01269  0.09988  0.19079 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -19.1642     7.4161  -2.584 0.013422 *  
## lnSuhu         8.6217     2.2311   3.864 0.000389 ***
## lnCurahHujan  -0.1766     0.1194  -1.479 0.146726    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1744 on 41 degrees of freedom
## Multiple R-squared:  0.3161, Adjusted R-squared:  0.2827 
## F-statistic: 9.473 on 2 and 41 DF,  p-value: 0.0004149

Model yang terbentuk dari summary di atas adalah Produktivitas = -19.1642 + 8.6217 Suhu - 0.1766 CurahHujan. Pada pemodelan diketahui nilai R-squared yang menunjukan 28.27% data yang dapat dijelaskan oleh model ini.

3.1 nilai ekspetasi

Eks <- predict(model_ekspektasi, newdata=data)

Eksp <- mean(Eks)
Eks
##        1        2        3        4        5        6        7        8 
## 8.151511 8.352887 8.288949 8.384564 8.195013 8.254872 8.249110 8.349758 
##        9       10       11       12       13       14       15       16 
## 8.314209 8.285576 8.350646 8.284466 8.320639 8.324290 8.280744 8.353113 
##       17       18       19       20       21       22       23       24 
## 8.327079 8.357225 8.522944 8.283143 8.313991 8.376012 8.445346 8.422967 
##       25       26       27       28       29       30       31       32 
## 8.415451 8.443594 8.411701 8.411047 8.519389 8.477118 8.589851 8.463893 
##       33       34       35       36       37       38       39       40 
## 8.438861 8.559024 8.499853 8.530978 8.582988 8.476459 8.506806 8.518851 
##       41       42       43       44 
## 8.495896 8.595122 8.648073 8.508991

3.2 Standar Deviasi

sd_eks <- sd(predict(model_ekspektasi, type = "response")) 
cat("Nilai standar deviasi dari G (std[G]) adalah",sd_eks, "\n")
## Nilai standar deviasi dari G (std[G]) adalah 0.1157838

Dari hasil tersebut diketahui bahwa nilai dari standar deviasi model ekspetasi Cobb Douglas adalah sebesar 0.115783

intercept <- model_ekspektasi$coefficients[1]
beta_suhu <- model_ekspektasi$coefficients[2]
beta_curahhujan <- model_ekspektasi$coefficients[3]

lnSuhu_range <- seq(min(data$lnSuhu), max(data$lnSuhu), length.out = 50)
lnCurahHujan_range <- seq(min(data$lnCurahHujan), max(data$lnCurahHujan), length.out = 50)

productivity_grid <- outer(lnSuhu_range, lnCurahHujan_range, function(x, y) {
  exp(intercept + beta_suhu * x + beta_curahhujan * y)
})

plot_ly(data, x = ~lnSuhu, y = ~lnCurahHujan, z = ~Produktivitas, type = "scatter3d", mode = "markers", marker = list(size = 3)) %>%
  add_surface(x = lnSuhu_range, y = lnCurahHujan_range, z = productivity_grid, colorscale = list(c(0, 1), c("blue", "green")), cmin = min(data$Produktivitas), cmax = max(data$Produktivitas)) %>%
  layout(title = "3D Scatter Plot ", scene = list(xaxis = list(title = "Log(Suhu)"), yaxis = list(title = "Log(Curah Hujan)"), zaxis = list(title = "Produktivitas")))
## Warning: 'surface' objects don't have these attributes: 'mode', 'marker'
## Valid attributes include:
## '_deprecated', 'autocolorscale', 'cauto', 'cmax', 'cmid', 'cmin', 'coloraxis', 'colorbar', 'colorscale', 'connectgaps', 'contours', 'customdata', 'customdatasrc', 'hidesurface', 'hoverinfo', 'hoverinfosrc', 'hoverlabel', 'hovertemplate', 'hovertemplatesrc', 'hovertext', 'hovertextsrc', 'ids', 'idssrc', 'legendgroup', 'legendgrouptitle', 'legendrank', 'lighting', 'lightposition', 'meta', 'metasrc', 'name', 'opacity', 'opacityscale', 'reversescale', 'scene', 'showlegend', 'showscale', 'stream', 'surfacecolor', 'surfacecolorsrc', 'text', 'textsrc', 'type', 'uid', 'uirevision', 'visible', 'x', 'xcalendar', 'xhoverformat', 'xsrc', 'y', 'ycalendar', 'yhoverformat', 'ysrc', 'z', 'zcalendar', 'zhoverformat', 'zsrc', 'key', 'set', 'frame', 'transforms', '_isNestedKey', '_isSimpleKey', '_isGraticule', '_bbox'

4 Model Variansi Cobb-Douglas

# Ekstraksi koefisien
beta_suhu <- model_ekspektasi$coefficients["lnSuhu"]
beta_curahhujan <- model_ekspektasi$coefficients["lnCurahHujan"]

# Menghitung variansi dan kovarians
var_suhu <- var(data$lnSuhu)
var_curahhujan <- var(data$lnCurahHujan)
cov_suhu_curahhujan <- cov(data$lnSuhu, data$lnCurahHujan)

# Menghitung kontribusi variansi
variance_contribution_suhu <- beta_suhu^2 * var_suhu
variance_contribution_curahhujan <- beta_curahhujan^2 * var_curahhujan
covariance_contribution <- 2 * beta_suhu * beta_curahhujan * cov_suhu_curahhujan

# Menghitung total variansi yang diverifikasi
total_variance_verified <- variance_contribution_suhu + variance_contribution_curahhujan + covariance_contribution

# Menampilkan hasil
cat("Kontribusi Variansi dari Suhu:", variance_contribution_suhu, "\n")
## Kontribusi Variansi dari Suhu: 0.01075106
cat("Kontribusi Variansi dari Curah Hujan:", variance_contribution_curahhujan, "\n")
## Kontribusi Variansi dari Curah Hujan: 0.001575346
cat("Kontribusi Kovarians antara Suhu dan Curah Hujan:", covariance_contribution, "\n")
## Kontribusi Kovarians antara Suhu dan Curah Hujan: 0.001079491
cat("Total Variansi yang Diverifikasi:", total_variance_verified, "\n")
## Total Variansi yang Diverifikasi: 0.01340589

Varian dari suhu memberikan kontribusi sebesar 0.01075106 dengan kontribusi yang paling tinggi. Kontribusi curah hujan adalah 0.001575346 dan Kontribusi Kovarians antara Suhu dan Curah Hujan adalah 0.001079491. Lalu saat semua kontribusi itu ditambah makan akan dihasilkan variansi yang terverifikasi sebesar 0.01340589.

5 Premi Asuransi

5.1 Premi dengan Loading Faktor

faktor_loading <- seq(0.01, 0.10, by = 0.01)
premi <- data.frame(faktor_loading)
premi$p_load <- NA
premi$p_load_unlog <- NA

for (i in 1:length(faktor_loading)) {
  premi$p_load[i] <- (1 + premi$faktor_loading[i]) * Eks
  premi$p_load_unlog[i] <- exp(premi$p_load[i])
}
## Warning in premi$p_load[i] <- (1 + premi$faktor_loading[i]) * Eks: number of
## items to replace is not a multiple of replacement length
## Warning in premi$p_load[i] <- (1 + premi$faktor_loading[i]) * Eks: number of
## items to replace is not a multiple of replacement length
## Warning in premi$p_load[i] <- (1 + premi$faktor_loading[i]) * Eks: number of
## items to replace is not a multiple of replacement length
## Warning in premi$p_load[i] <- (1 + premi$faktor_loading[i]) * Eks: number of
## items to replace is not a multiple of replacement length
## Warning in premi$p_load[i] <- (1 + premi$faktor_loading[i]) * Eks: number of
## items to replace is not a multiple of replacement length
## Warning in premi$p_load[i] <- (1 + premi$faktor_loading[i]) * Eks: number of
## items to replace is not a multiple of replacement length
## Warning in premi$p_load[i] <- (1 + premi$faktor_loading[i]) * Eks: number of
## items to replace is not a multiple of replacement length
## Warning in premi$p_load[i] <- (1 + premi$faktor_loading[i]) * Eks: number of
## items to replace is not a multiple of replacement length
## Warning in premi$p_load[i] <- (1 + premi$faktor_loading[i]) * Eks: number of
## items to replace is not a multiple of replacement length
## Warning in premi$p_load[i] <- (1 + premi$faktor_loading[i]) * Eks: number of
## items to replace is not a multiple of replacement length
data.table(premi)
##     faktor_loading   p_load p_load_unlog
##              <num>    <num>        <num>
##  1:           0.01 8.233027     3763.206
##  2:           0.02 8.314542     4082.814
##  3:           0.03 8.396057     4429.565
##  4:           0.04 8.477572     4805.767
##  5:           0.05 8.559087     5213.919
##  6:           0.06 8.640602     5656.735
##  7:           0.07 8.722117     6137.159
##  8:           0.08 8.803632     6658.386
##  9:           0.09 8.885147     7223.880
## 10:           0.10 8.966663     7837.401

5.2 Premi Standar Deviasi

premi_sd <- data.frame(faktor_loading)
premi_sd$p_load <- (sd_eks * premi_sd$faktor_loading) + Eksp
premi_sd$p_load_unlog <- exp(premi_sd$p_load)
data.table(premi_sd)
##     faktor_loading   p_load p_load_unlog
##              <num>    <num>        <num>
##  1:           0.01 8.407590     4480.947
##  2:           0.02 8.408748     4486.138
##  3:           0.03 8.409905     4491.335
##  4:           0.04 8.411063     4496.539
##  5:           0.05 8.412221     4501.748
##  6:           0.06 8.413379     4506.963
##  7:           0.07 8.414537     4512.185
##  8:           0.08 8.415695     4517.412
##  9:           0.09 8.416852     4522.645
## 10:           0.10 8.418010     4527.885
# Create plots
ggplot() +
  geom_line(data = premi, aes(x = faktor_loading, y = p_load_unlog), color = "blue", linetype = "solid") +
  geom_line(data = premi_sd, aes(x = faktor_loading, y = p_load_unlog), color = "red", linetype = "dashed") +
  labs(title = "Perbandingan Premi Asuransi dan Premi Standar Deviasi",
       x = "Faktor Loading",
       y = "Premi") +
  scale_x_continuous(breaks = seq(0.01, 0.10, by = 0.01)) +
  scale_y_continuous(labels = scales::comma) +
  theme_minimal()

model diatas ini menunjukan bahwa dengan menghitung ekspetasi premi membuat hasil premi lebih tinggi jika dibandingkan dengan menggunakan premi standar deviasi.